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Abstract 

This  report  summarizes  work  that  was  performed  in  fiscal  years  2002  through  2004 
on  or  related  to  the  “Facility  Unsteadiness  Research”  project  sponsored  by  the  Air 
Force  Office  of  Scientific  Research  (AFOSR).  Many  ground  test  facilities  experience 
some  form  of  flow  unsteadiness  at  various  points  in  their  operational  envelope.  This 
unsteadiness  can  at  times  cause  significant  problems  in  the  quality  of  data  obtained, 
the  structural  integrity  of  the  facility,  or  both.  Conventional  simulation  techniques, 
which  do  quite  well  for  steady-state  simulations,  have  been  found  to  be  deficient  in 
many  cases  when  applied  to  unsteady  flows.  The  current  project  attempts  to  address 
these  problems.  This  work  has  focused  on  three  areas:  building  up  the  computational 
infrastructure  of  a  production  CFD  code  to  be  able  to  tackle  large  scale  unsteady  prob¬ 
lems,  developing  a  rigorous  code  verification  capability  in  order  to  ensure  continued 
reliability  of  the  results  in  the  face  of  extensive  code  modification,  and  the  running 
of  test  cases  leading  toward  the  ultimate  objective  of  being  able  to  simulate  full  scale 
facilities.  The  Detached  Eddy  Simulation  model  and  the  hybrid  model  of  Nichols  and 
Nelson  were  found  to  offer  improved  results  over  conventional  turbulence  models  for 
unsteady  flows.  The  Method  of  Manufactured  Solutions  was  employed  to  verify  several 
of  the  numerical  methods  available  in  the  Wind-US  code.  The  algorithm  improvements 
to  the  Wind-US  code  have  now  made  it  possible  to  begin  to  model  large  ground  test 
facilities  and  resolve  medium  and  large-scale  instabilities  within  them. 
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Chapter  1 

Introduction 


Currently,  state  of  the  art  production  CFD  codes  solve  the  Reynolds  Averaged  Navier- 
Stokes  (RANS)  equations  using  implicit  schemes  which  are,  typically,  second  order 
(or  sometimes  third)  in  space  and  first  order  (sometimes  second)  in  time.  Turbulent 
closure  is  typically  obtained  by  using  a  one  or  two-equation  model.  This  combination 
is  generally  quite  successful  for  steady  state  simulations.  As  more  unsteady  simulations 
have  been  attempted,  however,  the  limitations  of  the  these  codes  have  become  apparent. 
The  sources  of  difficulty  can  be  broken  into  two  basic  categories:  the  numerical  scheme 
employed,  and  the  turbulence  model.  In  order  to  be  successful,  the  current  project  must 
address  both  issues. 

In  the  area  of  of  numerical  scheme  problems,  the  major  limitation  lies  in  the  ten¬ 
dency  of  the  methods  used  in  production  solvers  to  excessively  dissipate  the  unsteady 
fluctuations  in  the  flow.  Furthermore,  even  when  dissipation  is  within  limits,  dispersion 
errors  may  be  significant,  resulting  in  improper  prediction  of  wave  shapes  and  prop¬ 
agation  speeds.  Thus,  the  challenge  is  to  find  a  scheme  which  is  stable  enough  to  be 
used  in  a  production  setting,  but  accurate  enough  to  correctly  propagate  the  unsteady 
disturbances  from  their  source  to  the  region  of  interest,  which  in  a  large  facility  may  be 
hundreds  of  feet  removed. 

Turbulence  model  limitations  typically  arise  because  most  production  model  were 
designed  and  calibrated  for  RANS  solutions.  Thus,  by  design,  these  schemes  attempt  to 
model  all  turbulent  fluctuations,  whatever  the  length  or  time  scale.  Obviously,  in  an  un¬ 
steady  simulation  of  a  large  scale  test  facility,  some  of  the  fluctuations  which  one  would 
like  to  resolve  are  turbulent.  Therefore,  RANS  models  tend  to  be  overly  dissipative  in 
these  situations,  resulting  in  an  inability  to  accurately  propagate  the  fluctuations  in  the 
flow.  In  recent  years,  several  researchers  have  attempted  to  overcome  this  limitation 
by  combining  RANS  turbulence  models  with  Large  Eddy  Simulation  (LES)  models, 
using  the  former  in  regions  where  grid  resolution  is  insufficient  for  capturing  the  un¬ 
steady  fluctuations  and  the  latter  in  denser  portions  of  the  flow  domain.  These  models 
are  finding  increased  application,  but  because  their  strengths,  weaknesses,  and  failure 
modes  have  not  been  fully  mapped,  they  are  not  yet  generally  used  in  production  CFD 
environments. 

The  aforementioned  shortcomings  can  be  a  substantial  barrier  to  the  successful 


1 


AEDC-TR-05-1 


simulation  of  unsteady  flows  in  ground  test  facilities.  For  example,  the  16-S  supersonic 
wind  tunnel  at  AEDC  has  a  known  problem  with  oscillating  flow  angularity  in  the  test 
section. 

There  is  an  ongoing  effort  to  re-activate  and  upgrade  the  tunnel.  Included  in  this 
effort  are  fixes  addressing  the  flow  quality  problem.  This  required  an  understanding  of 
the  physics  involved  in  the  observed  unsteadiness  so  as  to  be  able  to  prevent/correct  it 
in  a  reasonably  efficient  manner. 

Because  the  capability  to  simulate  the  large-scale  unsteadiness  in  a  facility  as  siz¬ 
able  as  the  16-S  tunnel,  the  investigators  working  on  the  flow  quality  problem  were 
forced  to  construct  a  subscale  model  of  portions  of  the  larger  tunnel.  While  they  believe 
that  they  have  correctly  identified  the  problem  and  have  a  solution  for  it,  the  process 
has  been  slow  and  expensive.  It  is  believed  that  a  capability  to  numerically  simulate 
the  facility  could  have  saved  much  time  and  money,  had  it  been  available.  The  cur¬ 
rent  work  is  an  attempt  to  develop  such  a  capability.  The  work  encompassed  several 
different  (albeit  related)  areas:  investigating  two  approaches  to  hybrid  RANS/LES  tur¬ 
bulence  modeling,  building  up  of  code  infrastructure  to  better  handle  unsteady  flows, 
developing  code  verification  capabilities,  and  running  test  cases  to  establish  the  ability 
of  the  code  to  capture  the  necessary  physics.  One  of  the  decisions  made  early  in  this 
work  was  that  a  great  deal  of  emphasis  was  to  be  placed  on  the  transfer  of  the  result¬ 
ing  technology  into  a  production  environment.  For  this  reason,  as  well  as  the  overall 
complexity  of  the  problem  being  tackled,  most  of  the  work  was  performed  using  a 
modified  production  flow  solver,  rather  than  the  more  conventional  (and  simpler,  if  less 
flexible)  research  codes.  Also,  much  of  this  work  has  been  funded  by  sources  other 
than  AFOSR,  but  it  is  discussed  here  in  relation  to  the  Facility  Unsteadiness  Research 
project. 
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Chapter  2 

Theory  of  Hybrid  RANS/LES 
Turbulence  Models 


The  first  major  thrust  of  the  Facility  Unsteadiness  research  project  was  to  investigate 
two  different  approaches  to  hybrid  RANS/LES  turbulence  modeling.  The  purpose  was 
to  begin  to  map  the  strengths  and  weaknesses  of  such  methods,  as  well  as  gaining 
experience  in  running  the  particular  models  chosen. 

Hybrid  RANS/LES  turbulence  models  can  be  approached  theoretically  in  many 
ways.  One  way  is  to  see  them  as  derivatives  of  the  multi-scale  turbulence  models 
developed  by  Hanjalic,  Launder,  and  Schiestel[l],  Kim  and  Chen[2],  and  Duncan  and 
Liou[3].  These  turbulence  models  all  have  the  form 


dkp  dkp  3  / /  Vf  \  dkp\ 

dt  7  dxj  dxj  \  \  J  dxj  J 

dt  ox  j  dxj  \\  aEpJdxjJ 

dk,  3^___3_/V  Vt  \  dkt  \ 

dt  7  dxj  dxj  \  \  Gict  J  dxj  J 

3e,  d^r  _  d  ( (  vt  \  det  \ 

dt  7  dxj  dxj  \  \  0£,  J  dxj  J 
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where 


Vr  =c. 


( kP  +  kt)" 


(2.5) 


Equations  2.1  and  2.2  correspond  to  the  large-scale  turbulence.  Equations  2.3  and 
2.4  represent  the  smaller  more  dissipative  turbulent  scales.  The  large-scale  dissipation, 
Ep,  is  generally  the  same  size  as  the  small-scale  dissipation,  er,  for  these  multi-scale 
models.  The  large-scale  turbulent  kinetic  energy,  kp,  is  generally  an  order  of  magnitude 
larger  than  the  small-scale  turbulent  scale  energy,  kt. 

These  equations  suggest  two  possibilities  for  extension  to  unsteady  flows.  A  first 
approach  would  be  to  solve  equations  2.3  and  2.4  for  the  small-scale  turbulent  kinetic 
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energy  ( k To  implement  this  method,  the  source  terms  in  equations  2.3  and  2.4  would 
require  a  blending  function  based  on  the  filter  function  to  allow  them  to  be  used  when 
the  turbulent  length  scale  is  below  the  grid  resolution.  This  approach  would  also  require 
some  formulation  for  the  large-scale  dissipation  e;).  This  approach  could  require  a  good 
deal  of  tuning  to  be  applicable  over  a  wide  range  of  problems. 

A  second  approach  would  be  to  solve  equations  2.1  and  2.2  to  find  the  large-scale 
turbulence  parameters.  The  large-scale  turbulent  kinetic  energy  ( kp )  could  then  be  spa¬ 
tially  filtered  assuming  that  a  relationship  between  turbulent  kinetic  energy  and  the 
turbulent  length  scale  could  be  established.  The  eddy  viscosity  for  the  mean  flow  equa¬ 
tions  would  then  be  calculated  using  Equation  2.5  with  the  filtered  turbulent  kinetic 
energy.  This  approach  requires  that  the  turbulence  equations  be  solved  with  the  unfil¬ 
tered  eddy  viscosity  and  the  mean  flow  equations  be  solved  with  the  filtered  value  of 
eddy  viscosity. 

Both  approaches  have  been  applied  to  unsteady  flows.  One  example  of  the  first 
form  is  a  Detached  Eddy  Simulation  (DES)  turbulence  model  based  on  the  Spalart- 
Allmaras  one-equation  turbulence  model[4]  proposed  by  Spalart  etal.[ 5],  The  standard 
Spalart-Allmaras  turbulence  model  contains  a  destruction  term  for  eddy  viscosity  that 
is  inversely  proportional  to  the  square  of  the  distance  from  the  wall  ( d ).  The  hybrid 
version  replaces  the  wall  distance  with 

d  =  min  (d,  Cdes A)  (2.6) 

where  Cdes  is  0.65  and  A  is  a  characteristic  grid  length  scale.  The  modified  destruction 
term  has  the  effect  of  decreasing  the  eddy  viscosity  in  regions  not  in  the  immediate 
vicinity  of  walls.  A  similar  modification  to  Menter’s  Shear  Stress  Transport  Model 
(SST)  two-equation  model[6]  was  proposed  by  Strelets[7]. 

The  model  developed  by  Nichols  and  Nelson[8]  uses  the  second  approach  discussed 
above.  Equations  are  solved  for  the  large  scale  turbulent  parameters,  and  then  an  es¬ 
timate  is  made  as  to  how  much  of  the  turbulence  is  being  resolved  and  how  much 
must  be  modeled.  In  practice,  the  hybrid  model  is  implemented  by  modifying  the  eddy 
viscosity: 

V;  AM tRANS  (1  ~  (2-7) 

In  the  above  equation,  MtRANS  is  the  standard  eddy  viscosity  from  the  base  RANS  tur¬ 
bulence  model.  For  the  current  work,  both  the  Menter  SST  model[9]  (RANS  regime 
cases)  and  the  Jones-Launder  k  —  e  model[10]  (LES  regime  cases)  have  been  adapted. 
The  other  eddy  viscosity  (MtLES)  is  computed  based  on  a  form  taken  from  an  LES  Ad¬ 
equation  model[l  1]: 

V,LES  =  min  (c^Ay/k^s,  MtRANS )  (2.8) 

For  the  current  work,  a  value  of  0.0854  was  used  for  the  LES  model  coefficient  (CfjLES). 
The  subgrid  turbulent  kinetic  energy  is  defined  as: 

k^LES  =  fdkRANS  (2.9) 

where  the  clipping  function  is  defined  as  a  function  of  the  predicted  turbulent  length 
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scale  and  the  local  grid  size: 


(2.10) 


It  should  be  noted  that  other  implementations  of  this  model  (notably  that  in  NX  AIR) 
have  the  local  grid  spacing  (A)  rather  than  the  local  Nyquist  wavelength  (2A)  in  the 
denomenator  of  the  above  equation.  This  has  the  effect  of  increasing  the  effect  of  the 
LES  component  and  decreasing  the  RANS  component  of  the  turbulent  viscosity.  The 
turbulent  length  scale  used  in  this  effort  is  defined  by 


lt  =  max 


k- 

krans 

Zrans 


(2.11) 


where  CO  is  the  local  mean  flow  vorticity.  This  length  scale  is  a  mixture  of  the  traditional 
turbulent  scale  definition  for  two-equation  turbulence  models  {kI{ANS/ Erans)  and  the 
definition  usually  associated  with  algebraic  turbulence  models  ( ^/v, /(>) ).  The  factor  of 
six  was  chosen,  such  that  the  two  components  of  the  length  scale  were  approximately 
equal  in  a  simple  test  case.  This  turbulent  length  scale  definition  could  be  easily  adapted 
to  other  types  of  turbulence  models.  The  cut-over  function  (A)  is  defined  as: 

A=i(l+tanh(2n(7rf-£)^  (2.12) 

In  practical  terms,  then,  the  approach  that  Nichols  and  Nelson  have  taken  in  devel¬ 
oping  the  hybrid  model  for  unsteady  turbulent  flows  can  be  described  as  filtering  the 
RANS  turbulence  models  such  that  the  eddy  viscosity  does  not  include  the  energy  of 
grid-resolved  turbulent  scales.  Therefore,  the  filtered  RANS  turbulence  model  may  be 
thought  of  as  a  subgrid  model  for  very  large  turbulent  eddies.  The  objective  in  choosing 
this  filter  is  that  it  should  not  degrade  the  performance  of  the  turbulence  model  when 
the  largest  turbulent  scales  present  are  below  the  resolution  of  the  grid,  as  is  gener¬ 
ally  the  case  in  current  aircraft  CFD  applications-certainly  the  case  near  viscous  walls 
and  often  so  elsewhere.  This  approach  should  be  viable  for  the  current  class  of  CFD 
flow  solvers  and  would  not  require  any  more  grid  resolution  than  the  time  averaged 
CFD  grids  used  currently  except  in  areas  where  finer  details  of  the  instantaneous  flow 
structures  are  desired. 

As  a  final  note  on  hybrid  model  theory,  it  should  be  pointed  out  that  others  take 
widely  different  views  of  hybrid  models  and  the  theory,  if  any,  behind  them.  Durbin[12], 
for  example,  sees  problems  with  the  DES  model  (and  others  like  it)  because  “there  is 
an  influence  of  the  ensemble  on  the  instantaneous,  which  could  not  really  happen."  At 
the  same  conference  where  Durbin’s  aforementioned  ideas  were  presented,  the  DES 
model  was  described  by  Piomelli  and  Gatski  (in  the  question-and-answer  period  fol¬ 
lowing  a  presentation  by  Piomelli[13])  as  being  an  LES  model  that  has  a  simple  wall 
model  based  on  a  RANS  turbulence  model.  Thus,  while  there  may  be  (and  the  current 
work  has  found)  practical  benefits  to  the  hybrid  approach,  as  of  this  writing  there  is  no 
broad  consensus  as  to  its  theoretical  basis. 
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Chapter  3 

Evaluation  of  Hybrid 
RANS/LES  Turbulence  Models 
Using  an  LES  Code 


3.1  Introduction 

The  numerical  schemes  commonly  used  in  production  flow  solvers  are  relatively  dis¬ 
sipative,  and  many  such  schemes  have  themselves  been  found  to  act  much  like  a 
turbulence  model  (this  observation  being  the  basis  of  “monotone  implicit  LES”-aka 
“MILES”).  This  complicates  the  task  of  analyzing  the  effect  of  a  given  turbulence 
model  since  it  is  not  always  a  straightforward  process  to  determine  how  much  of  what 
is  observed  comes  from  the  numerical  scheme  and  how  much  comes  from  the  turbu¬ 
lence  model.  The  purpose  of  the  work  presented  in  this  chapter  is  to  investigate  hybrid 
turbulence  models,  comparing  them  to  both  conventional  RANS  models  and  LES  mod¬ 
els  in  the  context  of  a  code  designed  for  LES.  The  numerical  scheme  used  in  this  code 
has  comparatively  little  dissipation  (relative  to  production  RANS  codes);  therefore,  the 
effects  of  the  various  turbulence  models  should  be  more  clearly  shown. 

The  immediate  benefit  to  the  Facility  Unsteadiness  Research  project  is  a  clearer  un¬ 
derstanding  of  the  strength’s  and  weaknesses  of  the  various  models  when  operating  in 
LES  mode.  If  this  segment  of  the  operating  envelope  is  not  well  handled,  the  prospects 
for  successful  simulations  of  large  scale  facilities  would  be  relatively  dim.  Fortunately, 
as  will  be  shown,  these  models  have  been  found  to  perform  surprisingly  well.  Perfor¬ 
mance  when  the  configuration  is  more  in  the  RANS  range  will  be  discussed  in  a  later 
chapter. 

3.2  Description  of  the  Turbulence  Models 

The  two  hybrid  turbulence  models  discussed  in  the  previous  chapter  have  been  inves¬ 
tigated:  a  hybrid  k  —  £  model  using  the  approach  from  Nichols  and  Nelson  [14]  and 
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Table  3.1:  Flow  Parameters  for  Spatial  Mixing-Layer  Case 


7o,  K 

,  kPa 

M  i 

m2 

Mc 

U\,  m/s 

v2 

U, 

P2 

Pi 

Si,  mm 

291.0 

314.0 

1.80 

0.51 

0.52 

479.5 

0.355 

0.638 

8.0 

the  DES  model  of  Spalart[5],  For  comparison,  results  are  presented  using  the  RANS 
models  upon  which  the  two  hybrid  models  are  based:  the  Jones-Launder  k  —  £  model 
[10]  (with  compressibility  correction  [15])  and  the  Spalart-Allmaras  model  [4],  Fi¬ 
nally,  three  LES  models  have  been  run:  the  Smagorinsky  model  (as  implemented  by 
Hixon[16])  and  two  versions  of  a  /.-equation  model-one  with  constant  coefficients  and 
another,  dynamic  coefficient  version!  11],  An  attempt  was  made  to  run  without  any 
explicit  turbulence  model  (as  in  MILES),  but  the  code  was  not  stable  under  these  con¬ 
ditions,  and  the  attempt  was  unsuccessful.  This  is  not  surprising,  since  the  scheme  is 
not  monotone  and  is  therefore  inappropriate  for  a  MILES  approach. 

3.3  Description  of  the  Test  Case 

The  test  case  used  to  evaluate  the  models  was  based  on  the  mixing-layer  experiments 
of  Samimy  and  Elliot[17]:  a  two-stream  planar  mixing  layer,  with  one  stream  super¬ 
sonic  and  the  other  subsonic.  Experimental  data  were  available  at  various  locations 
between  the  trailing  edge  of  the  splitter  plate  and  210  mm  downstream  from  there.  The 
domain  of  the  simulations  included  a  short  length  of  splitter  plate  and  extended  over 
600  mm  downstream.  Only  the  middle  75  mm  of  the  experimental  rig  (out  of  150 
mm)  was  simulated.  Since  the  side  walls  were  not  being  modeled,  periodic  boundary 
conditions  were  employed  in  the  spanwise  direction.  This  approach  had  two  princi¬ 
pal  advantages.  First,  it  removed  any  need  for  injecting  artificial  turbulence  at  these 
boundaries  (when  flow  entered  the  domain  through  them).  Second,  it  permitted  the  use 
of  spanwise  averaging  to  decrease  the  simulation  time  required  to  achieve  statistical 
steady  state  (estimated  at  five  to  ten  times  faster  than  without  averaging). 

The  flow  parameters  for  this  case  are  given  in  Table  3.1.  In  the  table,  the  subscript 
numbers  1  and  2  refer  to  the  high-  and  low-speed  streams,  respectively,  and  8i  refers 
to  the  boundary-layer  thickness  at  the  end  of  the  splitter  plate  on  the  high-speed  side. 
Results  from  three  grid  resolutions  are  presented  here:  81  x61  x  34,  141  x91  x50,  and 
181  X  121  x  66.  At  the  finest  grid  resolution,  this  provided  about  thirty  points  across 
the  mixing  layer  at  the  210  mm  station.  A  pseudo-turbulent  fluctuation  field  was  added 
at  the  incoming  boundaries  to  better  match  the  experimental  conditions. The  solutions 
were  run  for  30,000  time  steps  to  wash  initial  transients  out  of  the  region  of  interest. 
Then  statistics  were  taken  over  60,000  further  steps.  As  mentioned  above,  since  trans¬ 
verse  boundary  conditions  were  periodic  (making  this  a  homogeneous  flow  direction), 
spatial  averaging  was  used  to  achieve  statistically  steady-state  results  relatively  quickly. 
The  time  step  for  all  cases  was  5.0  x  10~8  seconds. 

The  code  employed  for  this  work  uses  an  explicit  MacCormack-type  predictor- 
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Table  3.2:  Mixing  Layer  Growth  Rate  Predicted  by  Various  Turbulence  Models 


HM 

Exp. 

Smag. 

&-eqn. 

Dyn.  k 

Hybrid 

DES 

k  —  E 

Spalart 

db/dx 

0.0742 

0.0802 

0.0547 

0.0576 

0.0578 

0.0610 

0.0585 

0.0318 

0.0476 

1 

1.080 

0.737 

0.776 

0.779 

0.822 

0.788 

0.429 

0.642 

corrector  scheme  that  is  fourth  order  in  space  (on  Cartesian  grids)  and  second  order 
in  time.  Additional  details  on  the  geometry,  statistical  methods  (including  sampling 
procedures),  and  numerical  methods  may  be  found  in  Nelson[18].  Note,  however, 
that  the  domain  of  the  previous  simulations  did  not  extend  as  far  either  upstream  or 
downstream  as  does  this  study.  The  extended  computational  domain  was  included 
because  studies  of  the  effect  of  boundary  location  (not  presented  here)  for  the  current 
work  indicated  that  it  improved  results  somewhat.  Also  note  that  the  current  work  uses 
the  conventional  Favre-averaged  LES  governing  equations,  while  the  previous  work 
did  not.  In  other  respects,  the  techniques  used  are  much  the  same. 


3.4  Results 


The  results  presented  here  are  split  into  three  categories.  First,  a  general  survey  of 
the  behavior  of  the  various  models  is  presented.  Next,  a  study  of  the  effects  of  grid 
resolution  on  the  results  from  the  hybrid  models  is  presented  (this  doubles  as  a  grid 
convergence  study).  Finally,  the  self-similarity  predicted  by  the  hybrid  models  is  ex¬ 
amined. 

The  variables  shown  in  the  plots  are  now  defined.  The  non-dimensional  streamwise 
velocity  is  defined  as  follows: 


U  —  U2 


U1-U2 


(3.1) 


where  77  represents  the  time-averaged  velocity  in  the  streamwise  direction  and  U \  and 
U2  are  the  freestream  velocities  of  the  incoming  high-  and  low-speed  streams,  respec¬ 
tively.  The  shear  layer  thickness  at  a  given  streamwise  station  is  then  defined  as  the 
distance  between  the  point  where  U*  =0.1  and  the  point  where  it  reaches  0.9.  The 
normalized  vertical  coordinate  variable  used  in  the  plots  is  computed  using  the  loca¬ 
tion  where  U*  =  0.5  as  the  origin  and  the  vorticity  thickness  as  the  normalizing  factor: 


§m  = 


ih-ih 

(fj) 

\  J  /  max 


The  streamwise  turbulent  intensity  is  defined  as: 


<7„  = 


(3.2) 


(3.3) 


and  the  lateral  turbulent  intensity  is  defined  as: 


v  — 


<7, 


(3.4) 
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The  mv  component  of  Reynolds  stress  is: 

R„v  =  uv  —  uv  (3.5) 

The  streamwise  velocity  skewness  is  computed  as  m,3/o3  and  the  flatness  as  m,4/cj4. 


3.4.1  Behavior  of  the  Various  Models 


In  this  section,  the  predictions  of  the  various  turbulence  models  are  compared  with 
each  other.  The  results  presented  here  are  taken  from  the  fine  grid  simulations  (with 
grid  dimensions  of  181  x  121  x  66). 

As  shown  in  Fig.  3.1,  all  the  models  predict  roughly  the  same  initial  growth  of  the 
mixing  layer.  As  one  reaches  the  first  experimental  data  point  (60  mm  downstream  of 
the  trailing  edge  of  the  splitter  plate),  in  some  cases  the  differences  between  models  ap¬ 
proach  12  percent  of  the  experimentally  determined  thickness.  At  the  last  experimental 
data  point  (210  mm  downstream),  the  maximum  difference  between  models  is  almost 
22  percent  of  the  experimental  value.  Given  the  lack  of  data  available  on  the  incom¬ 
ing  streams,  and  the  consequent  assumptions  made  in  the  simulations,  it  is  perhaps  not 
surprising  that  the  absolute  values  of  shear  layer  thickness  would  differ  from  the  exper¬ 
iment.  In  addition,  the  geometry  has  been  somewhat  simplified  (and  under-resolved  as 
well)  from  the  experimental  rig.  This  simplification  would  also  have  an  effect  on  the 
initial  portions  of  the  mixing  layer. 

Of  more  importance  is  the  predicted  shear  layer  growth  rate,  which  should  be  more 
or  less  independent  of  the  details  of  the  incoming  streams  once  the  mixing  layer  has 
achieved  self-similarity.  In  the  experiment,  the  flow  was  found  to  be  self-similar  by 
the  time  it  reached  150  mm  downstream  of  the  splitting  plate.  The  shear  layer  growth 
rate  has  been  computed  using  a  least-squares  curve  fit  of  the  data  between  the  150  and 
300  mm  stations.  The  results  are  shown  in  Table  3.2.  The  incompressible  growth  rate 
((^) .)  is  computed  using  the  form  suggested  by  Papamoschou  and  Roshko[19]: 


db\  _Qy~ru)(l+^) 

dx  )  j  l+ru^/s 


(3.6) 


Where  r\j  is  the  velocity  ratio  (U2/U 1),  s  is  the  density  ratio  (P2/P1),  and  the  value  for 
the  coefficient  (Co)  was  0.083,  as  suggested  by  Goebel  and  Dutton[20]  for  thickness 
measured  as  the  distance  between  the  two  points  where  the  normalized  streamwise 
velocity  is  0.1  and  0.9,  respectively.  Not  surprisingly,  the  RANS  models  {k  —  £  and 
Spalart)  perform  worst  of  all.  The  algebraic  LES  model  (Smagorinsky)  improves  on 
their  results  somewhat.  Providing  even  better  results  are  the  one-equation  LES  models 
(the  /.-equation  and  dynamic  /.'-equation  models)  and  the  DES  model,  all  of  which 
predict  roughly  the  same  growth  rate.  Finally,  the  best  prediction  is  obtained  using 
the  hybrid  model.  All  of  the  simulated  mixing  layers  had  a  significantly  lower  growth 
rate  than  the  experiment.  The  experimental  value  was  obtained  from  only  three  points, 
however,  which  might  lead  to  some  error.  It  should  also  be  noted  that  there  is  a  great 
deal  of  scatter  in  experimental  measurements  by  different  researchers.  For  convective 
Mach  numbers  ranging  from  0.4  to  0.6,  the  results  (from  multiple  sources)  presented  in 
Elliot[21]  show  experimentally  measured  growth  rates  (normalized)  ranging  from  less 


9 


AEDC-TR-05-1 


than  0.5  to  more  than  1.0.  Thus,  in  addition  to  whatever  numerical  problems  there  may 
be  with  the  simulations,  one  would  expect  a  relatively  large  experimental  uncertainty 
as  well. 

Normalized  streamwise  velocity  profiles  are  shown  in  Fig.  3.2.  These  data,  and 
all  other  transverse  profiles  (unless  otherwise  noted),  are  taken  from  a  station  210  mm 
downstream  of  the  splitter  plate.  All  of  the  models  do  a  good  job  of  capturing  the 
experimentally  observed  behavior.  On  the  low-speed  side  of  the  layer,  the  DES  and 
Smagorinsky  models  appear  to  deviate  somewhat  more  than  the  other  models  (though 
the  errors  are  in  opposite  directions).  On  the  high-speed  side,  the  RANS  k  —  e  model 
underpredicts  the  velocity  slightly,  relative  to  the  other  models. 

The  inability  of  the  RANS  turbulence  models  to  properly  capture  the  unsteady  fea¬ 
tures  of  this  mixing  layer  become  obvious  in  Fig.  3.3,  which  shows  the  predicted 
streamwise  turbulence  intensity  across  the  mixing  layer.  Neither  of  the  RANS  models 
captures  the  correct  intensity  levels,  and  a  qualitative  comparison  shows  that  the  shapes 
of  the  curves  are  distorted  as  well,  particularly  by  the  k  —  £  model.  In  contrast,  the  other 
models  capture  the  flow  features  quite  well.  Note  that  away  from  the  core  of  the  mix¬ 
ing  layer,  the  grid  resolution  is  insufficient  to  maintain  the  freestream  turbulence  levels, 
resulting  in  a  drop-off  in  the  profiles  at  the  edge  of  the  mixing  layer,  relative  to  the  ex¬ 
perimental  results.  Much  the  same  behavior  is  observed  in  the  Reynolds  cross-stress 
term  ( RIIV ),  although  the  deviation  from  experiment  is  somewhat  greater  for  all  the 
models  (Fig.  3.4).  This  is  to  be  expected,  even  apart  from  any  error  in  the  simulations, 
since  the  experimental  error  for  cross  terms  is  higher  than  that  for  normal  terms.  This 
error  becomes  even  greater  as  one  moves  to  higher-order  turbulence  statistics  [22], 

The  streamwise  velocity  fluctuation  skewness  and  flatness  are  shown  in  Figs.  3.5 
and  3.6.  Again,  the  hybrid  and  LES  models  are  much  more  accurate  than  the  unmodi¬ 
fied  RANS  models.  For  these  higher  moments  of  turbulence,  however,  the  limitations 
of  the  grid  resolution  become  more  apparent.  Beyond  the  core  region  of  the  mixing 
layer,  the  effects  of  the  small-scale  turbulent  structures  become  more  important,  and 
the  lack  of  adequate  grid  resolution  for  these  small  flow  features  results  in  exagger¬ 
ated  peaks  (relative  to  the  experimental  data)  and,  in  the  case  of  the  flatness,  wild 
fluctuations.  Within  the  mixing  layer,  however,  the  flow  is  sufficiently  resolved  that 
the  experimentally  observed  features  are  well  captured  by  the  LES,  DES,  and  hybrid 
models. 

The  above  results  show  that,  for  this  case,  the  DES  and  hybrid  k— E  models  do  as 
well  or  better  than  the  pure  LES  models.  Relative  to  the  unmodified  RANS  models, 
the  DES  and  hybrid  models  represent  a  clear  improvement. 

3.4.2  Grid  Sensitivity 

There  are  two  parts  to  the  issue  of  grid  sensitivity.  The  first  is  whether  the  results 
computed  on  a  grid  are  of  sufficient  resolution  that  the  solution  truly  represents  the 
capabilities  of  the  turbulence  models.  The  second  is  whether  there  are  any  significant 
differences  between  the  responses  of  these  models  to  changing  grid  resolution.  To  an¬ 
swer  these  questions,  the  results  from  the  DES  and  hybrid  models  at  three  different  grid 
resolutions  are  presented.  One  complicating  factor  to  note  is  that,  unlike  a  traditional 
RANS  solution,  an  LES  solution  will  never  completely  “converge”  as  one  increases 
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grid  resolution  (until  it  becomes  a  DNS  solution).  This  is  because,  as  more  resolution 
is  added,  smaller  features  are  resolved,  resulting  in  additional  information  entering  the 
analysis.  Thus,  although  mention  is  made  here  of  “grid  convergence,”  the  meaning  of 
such  a  thing  is  somewhat  ambiguous.  The  working  definition  employed  here  is  that  a 
solution  is  considered  “grid  converged”  if,  in  the  core  of  the  mixing  layer,  the  differ¬ 
ences  between  the  fine  and  medium  grid  results  are  on  the  order  of  (or  better  than)  the 
differences  between  experimental  results  at  neighboring  stations  (when  plotted  using 
similarity  variables).  In  addition,  one  expects  that  the  differences  between  the  fine  and 
medium  grids  should,  in  general,  be  smaller  than  those  between  the  medium  and  coarse 
grids. 

In  Fig.  3.7  one  can  see  that  grid  convergence  appears  to  have  been  achieved  for  the 
shear  layer  thickness  on  both  the  fine  and  medium  resolution  grids.  The  same  appears 
to  be  true  of  the  streamwise  velocity  profiles  (Fig.  3.8).  Even  the  coarsest  grid  appears 
to  capture  these  basic  features  of  the  flow  quite  well. 

The  second-order  turbulence  moments  show  more  variability  with  grid  resolution 
(Figs.  3.9,  3.10,  and  3.11).  Even  here,  however,  the  coarsest  grid  does  a  reasonable  job 
of  capturing  both  the  qualitative  trends  and  even  the  fluctuation  magnitudes.  The  effect 
of  resolution  is  somewhat  different  for  the  different  models  in  that  the  largest  changes 
typically  are  at  different  places  on  the  profile.  The  DES  model  is  slightly  less  sensitive 
to  grid  resolution  than  the  hybrid  model  (this  is  particularly  obvious  in  Fig.  3.10).  The 
hybrid  model,  however,  appears  to  be  slightly  more  accurate  on  the  fine  grid,  especially 
at  the  edges  of  the  mixing  layer. 

Much  the  same  behavior  is  observed  in  the  higher  moments  (Figs.  3.12  and  3.13), 
although  there  appears  to  be  no  clear  grid  convergence  advantage  for  the  DES  model. 
The  hybrid  model  still  arguably  maintains  an  edge  in  accuracy  at  the  fringes  of  the 
mixing  layer,  particularly  on  the  low-speed  side. 

It  thus  seems  safe  to  say  that  the  fine  grid  resolution  is  sufficient  to  provide  mean¬ 
ingful  results.  It  would  also  appear  that,  for  the  most  part,  the  DES  and  hybrid  models 
provide  comparable  capability  at  the  various  grid  resolutions. 

3.4.3  Predicted  Self-Similarity 

One  of  the  most  interesting  features  of  mixing  layers  is  their  self-similar  nature.  The 
question  arises,  therefore,  do  the  DES  and  hybrid  models  predict  self-similar  behavior? 
If  so,  does  one  do  any  better  than  the  other?  To  examine  this  issue,  the  results  from 
several  downstream  stations  (where  experimental  data  were  available  for  comparison) 
are  presented. 

As  seen  in  the  experiments,  the  streamwise  velocity  achieves  self-similarity  very 
early  in  the  flow  (Fig.  3.14).  Both  the  hybrid  k  —  £  model  and  the  DES  model  capture 
this  behavior  well.  Although  it  may  require  magnification  to  see  it,  the  DES  model 
appears  to  predict  self-similarity  much  earlier  in  the  flow  than  either  the  hybrid  model 
or  the  experiment;  even  the  data  at  the  60  mm  station  show  this  trait.  The  hybrid  model 
does  somewhat  better,  although  the  deviation  from  similarity  at  60  mm  is  less  than 
that  observed  in  the  wind  tunnel.  Of  course,  the  closer  one  gets  to  the  splitter  plate, 
the  more  the  solution  depends  on  how  the  incoming  boundary  layers  and  freestream 
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turbulence  are  handled.  Therefore,  any  conclusions  drawn  from  these  observations  are 
tentative  at  best. 

In  the  experiments,  Elliott  observed  that,  while  streamwise  velocity  showed  self¬ 
similarity  by  the  120  mm  station,  the  self-similarity  of  higher  moments  of  turbulence 
was  delayed  until  the  150  mm  station.  Much  the  same  behavior  can  be  observed  in  the 
computational  results,  as  shown  in  Fig.  3.15  for  the  streamwise  turbulence  intensity. 
Both  turbulence  models  do  show  the  120  mm  station  as  being  much  closer  to  the  “final” 
self-similar  profile  than  was  observed  in  the  experiment.  In  general,  however,  self¬ 
similarity  is  achieved  that  is  comparable  to  the  experiment  not  only  in  peak  magnitude 
and  in  the  shape  of  the  curve,  but  also  in  the  variations  observed  between  stations. 

Figures  3.16  and  3.17  show,  respectively,  the  lateral  turbulence  intensity  and  the  uv 
component  of  the  Reynolds  stress.  For  these  and  following  plots,  only  data  from  the 
three  stations  furthest  downstream  have  been  plotted.  As  before,  the  models  behave 
slightly  differently,  with  DES  doing  better  in  some  places  and  the  hybrid  model  better 
in  others. 

Figures  3.18  and  3.19  show  third  and  fourth  moments  of  the  streamwise  velocity 
fluctuations.  Both  models  to  a  good  job  of  capturing  the  observed  self-similar  behavior 
within  the  core  of  the  layer.  At  the  edges  of  the  layer,  deviations  from  the  experimental 
data  become  significant,  although  the  curve  shapes  remain  consistent  (except  on  the 
high-speed  side  of  the  flatness). 

The  above  results  show  that  both  the  DES  and  hybrid  k  —  £  models  capture  the 
important  self-similarity  features  of  the  flow.  Neither  one,  however,  shows  any  clear 
advantage  over  the  other. 

3.5  Conclusions 

In  conclusion,  it  has  been  shown  that  both  the  DES  and  hybrid  k  —  £  models  provide,  for 
unsteady  “FES-like”  flows,  a  radical  improvement  over  the  conventional  RANS  models 
from  which  they  are  derived.  Results  using  these  models  compare  favorably  with  those 
obtained  using  conventional  FES  models.  Beyond  that,  it  is  extremely  difficult  to  draw 
any  conclusions.  Testing  has  shown  that  the  results  may  be  sensitive  to  a  wide  variety  of 
factors,  including  boundary  condition  treatment  at  the  inflow  and  outflow  boundaries, 
treatment  of  incoming  freestream  turbulence,  and  grid  resolution,  as  well  as  statistical 
sampling  time  and  methodology.  Each  of  these  factors  has  been  investigated,  and  an 
attempt  has  been  made  to  reduce  their  effect,  but  the  uncertainties  they  introduce  are 
on  the  order  of  the  observed  differences  between  the  DES  and  hybrid  k  —  £  models.  In 
addition,  while  the  numerical  scheme  used  is  much  less  dissipative  than  that  typically 
used  in  production  RANS  solvers,  it  is  not  without  its  limitations.  Finally,  some  exper¬ 
imental  error  exists  as  well.  Therefore,  given  the  relative  equivalence  of  the  results  and 
the  uncertainties  involved,  it  is  not  clear  how  much  of  the  deviation  from  experimen¬ 
tal  results  observed  in  the  simulations  is  due  to  turbulence  model  limitations  and  how 
much  is  due  to  errors  unrelated  to  the  turbulence  models  (not  to  mention  experimental 
uncertainty).  These  questions  are  addressed  in  the  following  chapter. 
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Figure  3.1:  Mixing  layer  thickness 
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Figure  3.2:  Non-dimensional  streamwise  velocity 
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Figure  3.3:  Streamwise  turbulence  intensity 
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Figure  3.4:  Reynolds  stress  profiles 
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Figure  3.5:  Streamwise  velocity  fluctuation  skewness 
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Figure  3.6:  Streamwise  velocity  fluctuation  flatness 
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(a)  Hybrid  k  —  e  model 


(b)  DES  model 


Figure  3.7:  Effect  of  grid  resolution  on  mixing  layer  thickness  results 
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Figure  3.8:  Effect  of  grid  resolution  on  non-dimensional  streamwise  velocity 
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Figure  3.9:  Effect  of  grid  resolution  on  streamwise  turbulence  intensity 
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Figure  3.10:  Effect  of  grid  resolution  on  lateral  turbulence  intensity  profiles 
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Figure  3.11:  Effect  of  grid  resolution  on  resolved  Reynolds  stress 
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Figure  3.12:  Effect  of  grid  resolution  on  streamwise  turbulence  fluctuation  skewness 
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Figure  3.13:  Effect  of  grid  resolution  on  streamwise  turbulence  fluctuation  flatness 
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Figure  3.14:  Self-similarity  of  streamwise  velocity 
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Figure  3.15:  Self-similarity  of  streamwise  turbulence  intensity 
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Figure  3.16:  Self-similarity  of  lateral  turbulence  intensity 
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Figure  3.17:  Self-similarity  of  Reynolds  stress 
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Figure  3.18:  Self-similarity  of  streamwise  velocity  fluctuation  skewness 
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Figure  3.19:  Self-similarity  of  streamwise  velocity  fluctuation  flatness 
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Chapter  4 

Investigation  into  the  Precision 
of  Unsteady  Simulations 

4.1  Introduction 

As  the  data  presented  in  the  previous  chapter  make  manifestly  clear  (see  also  Nel¬ 
son  and  Nichols[23]),  simulations  using  the  DES  model  often  yield  results  which  are 
very  similar  to  those  obtained  using  the  Nichols  and  Nelson  hybrid  approach.  Since 
other  sources  of  error  are  known  to  exist  in  these  simulations,  the  question  arises  as 
to  whether  or  not  the  observed  differences  due  to  model  choice  are  actually  significant 
relative  to  the  other  error  sources.  If  not,  then  it  becomes  vital  that  work  on  unsteady 
flow  simulations  (of  this  nature)  focus  more  on  the  dominant  error  sources  and  less  on 
turbulence  model  issues. 

To  answer  this  question,  a  series  of  simulations  was  performed,  using  the  same  LES 
code  employed  for  the  previous  chapter’s  work,  which  examined  a  variety  of  issues 
to  determine  their  relative  effect  on  the  results.  The  test  case  was  again  the  planar 
mixing  layer  (convective  Mach  number  of  0.51)  studied  experimentally  by  Samimy 
and  Elliott[17,  24,  25], 

The  potential  sources  of  error  in  a  CFD  simulation  of  any  type  are  nearly  limitless. 
Therefore,  of  necessity,  the  following  analysis  only  examines  a  subset  of  issues  that 
could  be  influencing  the  results.  The  areas  investigated  can  be  broken  into  two  broad 
themes:  boundary  condition  issues  and  statistical  issues. 

Since  any  solution  to  the  Navier-Stokes  equations  is  determined  by  the  boundary 
conditions  imposed  in  space  and  time  (in  conjunction  with  the  physical  geometry  of 
the  case),  the  treatment  of  boundaries  is  an  obvious  area  to  investigate  for  impact  on 
the  overall  solution.  This  is  particularly  true  for  unsteady  cases,  where  the  reflection 
of  transient  waves  from  boundaries  can  be  a  serious  problem.  Thus,  the  placement  of 
both  the  inflow  and  outflow  boundaries  was  varied  in  an  attempt  to  determine  the  extent 
to  which  this  influences  the  solution.  Also,  the  original  inflow  boundary  condition 
employed  for  these  simulations  was  based  on  an  inviscid  (characteristic)  analysis.  To 
help  assess  the  effect  of  such  an  assumption,  the  boundary  condition  was  modified  so 
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as  to  add  the  viscous  stresses  as  well. 

Another  possible  source  of  error  for  unsteady  simulations  is  the  procedure  by  which 
the  unsteady  solution  is  processed  to  obtain  time-averaged  results.  If  the  period  over 
which  the  solution  is  averaged  is  not  sufficient,  then  the  results  will  not  be  accurate, 
particularly  for  the  higher  statistical  moments.  Furthermore,  if  the  solution  has  not 
been  given  sufficient  time  to  wash  out  the  effects  of  the  initial  conditions,  then  these 
residual  effects  can  also  contaminate  the  solution.  The  current  analysis  attempts  to 
address  both  of  these  issues. 

Among  the  issues  not  addressed  in  this  work  were  the  effects  of  the  simplified  ge¬ 
ometry  used  for  these  simulations  and  the  specifics  of  the  incoming  pseudo-turbulence 
(spectrum,  amplitude,  etc.).  While  these  issues  might  be  important  at  times,  these  ar¬ 
eas  were  ignored  in  order  to  focus  on  the  aforementioned  areas.  This  was  partly  due  to 
the  project’s  limited  resources,  but  also  because  of  the  lack  of  experimental  data  with 
which  to  compare  the  simulated  inflows,  especially  on  the  low  speed  side. 

4.2  Test  Case  Description 

The  effect  of  grid  resolution  on  the  issues  under  consideration  was  tested  by  running 
the  tests  using  three  different  base  grids.  The  resolutions  for  these  base  grids  were 
91  x  61  x  34,  121  x  91  x  50,  and  161  x  121  x  66.  By  default,  the  same  boundary 
conditions  were  used  and  the  same  statistical  data  collection  procedure  was  followed 
(run  for  30,000  steps  and  then  collect  data  for  a  further  60,000).  Except  where  noted,  a 
variation  of  the  Nichols  and  Nelson  hybrid  model  was  used  for  all  tests. 

The  boundary  location  tests  varied  the  above  grids  as  follows  (see  Table  4.1).  The 
first  modification  was  to  extend  the  inflow  by  ten  uniform  points  in  the  upstream  di¬ 
rection.  Second,  for  the  extended  outflow  case,  a  region  of  ten  highly  stretched  points 
was  added  to  the  outflow  boundary.  Then  a  third  grid  was  constructed  with  both  of  the 
aforementioned  extensions. 

A  fourth  grid  was  constructed  similar  to  the  extended  inflow  case,  but  with  a  further 
twenty  points  added  to  the  inflow,  for  a  total  of  thirty  more  streamwise  points  than  the 
base  grid.  The  fifth  grid  was  similar  to  the  extended  outflow  grid,  but  with  an  additional 
ten  highly  stretched  points  at  the  outflow.  The  sixth  and  final  geometric  variation  had 
both  the  additional  inflow  and  outflow  points.  The  effects  of  adding  viscous  terms  to 
the  inflow  boundary  were  then  tested  using  the  “Ext.  in/out-2”  grid. 

To  test  the  precision  of  the  statistical  data  collection  procedure,  three  additional 
cases  were  run  (again  using  the  “Ext.  in/out-2”  grid).  Two  of  these  cases  extended  the 
length  of  time  over  which  data  was  collected.  The  first  of  these  took  data  for  90,000 
iterations,  rather  than  the  base  60,000;  the  second  case  doubled  the  base  period  for 
collecting  data  to  120,000  time  steps.  The  final  variation  delayed  taking  statistical  data 
until  the  solution  had  run  for  90,000  steps;  at  that  point  data  was  collected  for  the 
standard  60,000  steps. 
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Table  4. 1 :  Unsteady  Simulation  Precision  Test  Cases 


Label 

Description 

Base 

91  x  61  x  34, 141  x  91  x  50,  or  161  x  121  x  66. 

Extend  in 

10  uniform  points  added  upstream 

Extend  out 

10  stretched  points  added  downstream 

Extend  in/out 

Extend  both  inflow  and  outflow  10  points 

Ext.  in-2 

30  uniform  points  added  upstream 

Ext.  out-2 

20  stretched  points  added  downstream 

Ext.  in/out-2 

Add  30  pts  upstream  and  20  pts  downstream 

Vise,  inflow 

Add  viscous  stresses  at  inflow 

Ext.  stats 

Extend  time  averaging  by  50% 

Ext.  stats-2 

Take  data  for  twice  as  long  as  base 

Ext.  init. 

Delay  start  of  averaging 

4.3  Results 

As  was  seen  before  [23],  the  non-dimensional  streamwise  velocity  profile  was  com¬ 
pletely  insensitive  to  all  variations.  An  example  of  this  is  plotted  in  Figure  4.1,  which 
shows  the  profile  at  the  210mm  station  at  the  coarsest  grid  resolution  (base  grid  91  x 
61  x  34).  In  addition  to  the  aforementioned  test  cases,  the  previously  presented  results 
from  the  DES  model  are  plotted  for  comparison.  In  this  work,  the  choice  of  model 
is  considered  to  be  “significant”  only  if  the  DES  model  results  lie  outside  the  scatter 
from  the  other  cases.  In  this  case,  obviously,  the  results  are  insensitive  to  the  choice  of 
model. 

A  more  interesting  behavior  is  found  in  the  streamwise  turbulent  intensity.  As 
shown  in  Figure  4.2,  the  relative  importance  of  the  turbulence  model  varies  with  grid 
resolution.  At  the  coarsest  resolution,  the  DES  model  prediction  is  noticeably  different 
on  the  high  speed  side  of  the  mixing  layer.  At  the  medium  resolution,  however,  this 
difference  is  reduced  to  the  point  of  insignificance.  For  the  fine  grid  cases,  the  DES 
model  again  appears  to  result  in  significant  (albeit  small)  differences.  It  is  encouraging 
that,  for  all  three  baseline  grid  resolutions,  the  boundary  location  does  not  dramatically 
change  the  character  of  the  results.  Therefore,  while  it  may  be  difficult  to  determine 
which  of  the  hybrid  approaches  is  best,  the  dramatic  improvement  observed  relative  to 
the  base  RANS  models  (see  Figure  3.3)  appears  to  be  real. 

A  stronger  dependence  on  boundary  location  is  shown  by  the  higher  moments  of 
turbulence.  This  is  illustrated  in  Figure  4.3  using  the  streamwise  velocity  skewness. 
In  this  case,  choice  of  hybrid  model  seems  to  be  less  important,  with  the  scatter  from 
the  boundary  location  variations  bracketing  the  DES  curves  virtually  everywhere.  As 
expected,  the  largest  deviations  occur  at  the  edges  of  the  mixing  layer,  while  the  core 
shows  good  agreement  (improving  with  resolution)  for  all  grids.  Of  course,  given  the 
limited  grid  resolution  beyond  the  core  of  the  mixing  layer,  it  is  doubtful  that  any  strong 
conclusions  can  be  formed  based  on  that  deviation.  Again,  the  differences  observed  are 
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not  such  that  the  overall  benefit  of  hybrid  models  relative  to  the  baseline  RANS  models 
would  be  called  into  question. 

The  effects  of  the  variation  of  the  statistical  sampling  procedure  or  inclusion  of 
viscous  terms  at  the  inflow  boundary  are  shown  in  Figure  4.4  for  the  streamwise  turbu¬ 
lent  intensity.  Within  the  core  of  the  mixing  layer,  only  delaying  the  initiation  of  data 
collection  has  an  appreciable  effect  (with  a  reduction  in  peak  intensity  at  the  medium 
resolution).  This  result  was  unexpected,  especially  since  extending  the  statistical  gath¬ 
ering  period  (so  that  the  two  runs  were  of  the  same  length)  had  no  such  effect.  Neither 
the  fine  grid  nor  the  coarse  grid  showed  the  same  behavior.  A  similar  situation  was 
observed  in  the  uv  component  of  Reynolds  stress  (Figure  4.5).  In  general,  extending 
the  period  before  gathering  statistics  had  a  larger  than  expected  effect  on  the  low  speed 
side  of  the  mixing  layer,  while  the  core  and  high  speed  sides  were  relatively  unaffected. 
This  is  well  illustrated  by  the  streamwise  velocity  skewness,  as  shown  in  Figure  4.6. 

In  the  above  plots,  the  addition  of  viscous  terms  at  the  inflow  was  observed  to  have 
little  effect.  When  changes  were  observed,  they  were  confined,  for  the  most  part,  to 
the  edges  of  the  mixing  layer.  In  this  region,  the  other  uncertainties  in  the  simulation 
make  it  difficult  to  draw  any  definite  conclusions  as  to  advantages  or  disadvantages  of 
this  strategy. 

4.4  Conclusions 

Based  on  the  above  analysis,  it  appears  that  the  question  of  whether  or  not  the  choice 
of  hybrid  model  technique  is  significant  depends  to  some  degree  on  what  data  is  sought 
and  the  relative  resolution  of  the  grid.  While  the  precision  of  these  simulations  is 
clearly  sufficient  to  demonstrate  the  general  superiority  of  the  hybrid  technique  over 
a  standard  RANS  model  for  this  class  of  problem  (as  in  Nelson  and  Nichols[23]), 
the  uncertainties  in  these  results  would  appear  to  call  into  question  most  conclusions 
regarding  the  relative  merit  of  the  different  hybrid  techniques.  For  most  quantities,  the 
differences  are  either  not  significant  at  all,  or  if  they  are,  it  is  not  yet  clear  which  is  the 
better  model. 

Furthermore,  while  a  mixing  layer  is  expected  to  behave  in  relatively  “universal” 
fashion  once  self-similarity  is  achieved  (and  the  data  is  plotted  using  the  appropriate 
similarity  variables),  the  initial  development  is  heavily  influenced  by  the  inflow  condi¬ 
tions  (e.g.  boundary  layer  thickness,  inflow  turbulence  intensity,  turbulence  spectrum, 
etc.).  It  is,  therefore,  entirely  possible  (even  probable),  that  even  if  all  of  the  other  issues 
raised  here  were  dealt  with,  the  inflow  conditions  could  be  tweaked  to  give  an  advan¬ 
tage  to  either  of  the  hybrid  models.  In  the  absence  of  more  complete  experimental  data 
about  the  inflow  conditions,  however,  such  an  exercise  would  prove  nothing. 

As  with  boundary  location,  this  investigation  found  that  the  original  statistical  data 
collection  methodology  was  adequate  to  clearly  show  the  advantages  of  the  hybrid 
models  to  predict  the  core  of  the  mixing  layer  relative  to  the  RANS  models  upon  which 
they  are  based.  Extending  the  sampling  time  or  delaying  the  initiation  of  data  sampling 
was,  however,  seen  to  have  an  unexpectedly  large  effect.  Given  that  the  preponderance 
of  changes  were  seen  on  the  low-speed  side  of  the  mixing  layer,  it  would  appear  that 
the  original  sampling  technique  did  not  allow  for  the  core  flow  on  that  side  to  develop 
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Figure  4.1:  Variation  of  non-dimensional  streamwise  velocity  with  boundary  location 
(coarse  grid) 


sufficiently  to  ensure  a  statistically  stationary  solution  everywhere  in  the  flowfield. 

Adding  in  viscous  terms  at  the  inflow  boundary  was  seen  to  have  only  a  minor 
effect  in  most  areas  of  the  flow.  When  differences  were  observed,  it  was  not  readily 
apparent  that  the  results  were  thereby  improved.  It  is  expected  that  such  a  question  can 
only  be  resolved  using  simpler  test  cases  which  have  a  more  precisely  known  solution. 
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(a)  Coarse  base  grid 


(b)  Medium  base  grid 


(c)  Fine  base  grid 


Figure  4.2:  Effect  of  boundary  location  on  non-dimensional  streamwise  turbulent  in¬ 
tensity 
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Figure  4.3:  Effect  of  boundary  location  on  non-dimensional  streamwise  velocity  skew¬ 
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(a)  Coarse  base  grid 


(b)  Medium  base  grid 


(c)  Fine  base  grid 


Figure  4.4:  Variation  of  streamwise  turbulent  intensity  with  data  collection  method 
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Figure  4.5:  Variation  of  uv  component  of  Reynolds  stress  with  data  collection  method 
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Figure  4.6:  Variation  of  streamwise  velocity  skewness  with  data  collection  method 
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Chapter  5 

Applications  of  Hybrid 
RANS/LES  Turbulence  Models 


The  standard  Spalart-Allmaras  one-equation  turbulence  model  and  Menter’s  SST  two- 
equation  turbulence  model  along  with  the  Detached  Eddy  Simulation  (DES)  version  of 
the  Spalart-Allmaras  model  and  the  hybrid  version  of  the  SST  model  were  applied  to 
Several  test  cases.  The  first  two  cases  are  two-dimensional  flows.  In  reality  there  are  no 
truly  two-dimensional  unsteady  turbulent  flows,  but  these  two-dimensional  cases  allow 
the  turbulence  models  to  be  investigated  quickly  and  allow  some  means  of  characteriz¬ 
ing  their  performance.  Three-dimensional  effects  would  tend  to  reduce  the  magnitude 
of  the  unsteady  forces  and  delay  the  onset  of  the  unsteadiness  in  these  cases  since  the 
third  dimension  would  provide  an  avenue  to  relieve  the  maximum  and  minimum  pres¬ 
sure  peaks.  The  remaining  test  cases  were  true  three-dimensional  applications.  All 
these  test  cases  were  run  using  the  NXAIR  Navier-Stokes  code[26].  Wall  function 
boundary  conditions  were  used  for  these  applications  since  they  have  been  shown  to 
provide  reasonable  results  for  similar  unsteady  applications  [26]. 


5.1  Two-Dimensional  Cylinder 

Unsteady  two-dimensional  calculations  were  performed  for  the  vortex  shedding  from 
a  circular  cylinder  for  M  =  0.2  and  Re,]  =  8x  106.  The  computational  grid  was  201  x 
201  with  an  initial  wall  spacing  of  2  x  1 0  4  diameters,  which  corresponds  to  a  y+  of 
about  20.  The  time  step  used  in  the  calculations  was  a  constant  9.1  x  1 0  5  sec.,  which 
corresponds  to  about  150  steps  per  cycle  of  the  shedding  frequency. 

The  instantaneous  pressure  contours  for  the  four  turbulence  models  are  shown  in 
Fig.  5.1.  The  vortex  shedding  is  clearly  evident  behind  the  cylinder.  The  hybrid 
RANS/LES  models  produce  a  much  lower  pressure  in  the  core  of  the  shed  vortices 
than  do  the  standard  models.  The  pressure  distribution  across  the  nearest  shed  vortex 
for  each  turbulence  model  is  shown  in  Fig.  5.2  (note  that  the  vortex  location  is  not 
exactly  the  same  for  each  of  the  models  at  this  observation  time).  The  legend  notation 
is  as  follows:  “SST”  refers  to  the  standard  Menter  SST  model,  “SST  MS”  refers  to  the 
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hybrid  SST  model,  “SA”  refers  to  the  Spalart-Allmaras  model,  and  “SA  DES”  refers 
to  the  DES  version  of  the  Spalart-Allmaras  model.  The  hybrid  RANS/LES  models  are 
seen  to  predict  about  the  same  minimum  pressure  in  the  vortex.  This  would  indicate 
an  improvement  in  the  results  with  the  hybrid  RANS/LES  models  since  the  standard 
RANS  models  are  known  to  over-predict  the  minimum  pressure  in  a  vortex[8].  The 
instantaneous  eddy  viscosity  contours  are  shown  in  Fig.  5.3.  The  eddy  viscosity  distri¬ 
butions  differ  radically  for  each  turbulence  model.  The  SST  model  predicts  the  highest 
eddy  viscosity  and  also  has  the  highest  pressure  in  the  vortex  core.  The  Spalart  model 
has  high  eddy  viscosity  near  the  cylinder  and  rapidly  decreases  away  from  the  body. 
The  hybrid  SST  model  has  low  eddy  viscosity  near  the  cylinder,  but  the  eddy  viscosity 
rapidly  increases  away  from  the  cylinder  as  the  computational  grid  becomes  coarser. 
The  DES  model  has  relatively  low  eddy  viscosity  everywhere  in  the  field.  The  dif¬ 
ferences  in  the  DES  model  and  the  hybrid  SST  model  may  be  due  to  the  fact  that  the 
latter  examines  locally  predicted  turbulent  length  scales  when  determining  whether  to 
operate  in  RANS  or  LES  mode,  while  the  former  switches  based  purely  on  distance 
from  the  nearest  wall  and  the  local  grid  scale. 

The  spectral  results  for  axial  and  normal  force  are  shown  in  Fig.  5.4  and  Fig. 
5.5  respectively.  These  results  were  taken  from  the  last  4096  time  steps  of  the  solu¬ 
tion.  While  the  energy  in  the  primary  peak  is  similar  for  all  the  models,  the  hybrid 
RANS/LES  models  are  seen  to  produce  much  more  energy  in  the  secondary  peaks  and 
the  overall  background.  The  non-harmonic  higher  energy  peaks  at  the  higher  frequen¬ 
cies  predicted  by  the  hybrid  RANS/LES  models  indicate  the  presence  of  small-scale 
structures  in  the  flow  field.  These  structures  can  clearly  be  seen  in  Fig.  5.1.  The  av¬ 
erage  cylinder  drag  and  the  Strouhal  number  of  the  primary  peak  in  the  normal  force 
spectrum  for  both  the  CFD  and  experimental  data[27,  28,  29]  are  shown  in  Table  5.1. 
The  predicted  Strouhal  number  is  seen  to  be  in  good  agreement  with  the  experimental 
data  for  all  the  turbulence  models.  The  hybrid  RANS/LES  models  predict  a  higher 
drag  coefficient  than  do  the  standard  turbulence  models,  but  all  the  models  produce  a 
reasonable  result  compared  to  the  scatter  in  the  experimental  values. 

The  data  from  Jones  [27]  was  taken  in  a  16  foot  tunnel  with  a  3  foot  cylinder.  This 
represents  a  tunnel  blockage  of  nearly  20  percent.  Therefore,  the  effect  of  the  walls 
is  potentially  quite  large.  To  test  this,  the  case  was  re-run  with  walls  present  using 
the  standard  SST  model.  With  the  walls  in  place  the  results  from  this  run  indicate  a 
significant  effect  on  the  predicted  drag.  Indeed,  the  drag  now  matches  the  experimental 
data  (t'd  =  0.53). 

5.2  Two-dimensional  Oscillating  NACA  0015  Airfoil 

The  second  test  case  was  an  oscillating  two-dimensional  NACA  0015  airfoil.  The 
experimental  configuration  [30]  included  a  wing  section  that  spanned  the  wind  tunnel 
wall  and  an  end  plate.  Pressure  data  was  taken  at  four  spanwise  locations  and  did 
indicate  that  three-dimensional  effects  were  not  insignificant  at  the  post-stall  angles-of- 
attack.  This  could  also  be  seen  in  tuft  photographs.  Sectional  coefficients  were  taken 
by  integrating  the  pressure  measurements.  The  comparisons  here  are  with  the  mid-span 
location  that  had  twenty  pressure  ports.  The  pressure  transducer  output  was  low-pass 
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Table  5.1:  Comparison  of  CFD  and  Data  for  the  Unsteady  Flow  from  a  Circular  Cylin¬ 
der  for  M  =  0.2  and  Red  =  8.0  x  106 


Data  Source 

Cd 

St 

CFD-SST  Model 

0.62 

0.29 

CFD-Hybrid  SST  Model 

0.85 

0.30 

CFD-Spalart-Allmaras  Model 

0.68 

0.29 

CFD-DES  Model 

0.83 

0.30 

Data-Jones  et.  al  [27] 

0.53 

0.30 

Data-Roshko  [28] 

0.79 

0.27 

Data-Schlichting  [29] 

- 

0.29 

filtered  at  500  Hz.  The  net  effect  of  the  experimental  setup  was  that  the  unsteady  data 
is  highly  filtered. 

The  flow  conditions  for  the  experiment  was  M  =  0.29  and  Rec  =  1 .95  x  106,  where 
c  is  the  airfoil  chord.  The  airfoil  was  oscillated  about  the  quarter  chord  location  at  three 
mean  angles-of-attack,  4°,  11°,  and  15°.  The  amplitude  of  the  oscillation  was  ±4.2° 
for  all  cases.  The  instantaneous  angle-of-attack  was  given  by 

a(t)  =  ao  +  ai  sin  (Kt)  (5.1) 

where  K  is  the  reduced  frequency  given  by 

CDc 

K  = - =  2  St  (5.2) 

Here  CO  was  the  frequency  of  oscillation  (10.1  Hz)  which  yields  K  =  0.1. 

The  calculations  were  performed  on  a  341  x  51  "C"  mesh  with  241  points  around 
the  airfoil.  The  initial  wall  spacing  was  0.00025  chords  that  corresponds  to  a  y+  of 
about  30.  Wall  function  boundary  conditions  were  used  in  these  calculations.  The  time 
step  for  the  calculation  was  2.2  x  1 0  5  seconds  which  corresponds  to  about  3200  time 
steps  per  wing  oscillation  cycle.  The  pressures  around  the  wing  were  integrated  at  each 
time  step  for  comparison  with  the  experimental  data. 

Figure  5.6  shows  the  pressure  contours  for  the  airfoil  oscillating  about  the  4°  mean 
angle-of-attack.  The  airfoil  is  at  a  =  7.92°  on  the  down  stroke.  There  appears  to  be 
very  little  difference  between  the  models  for  this  attached  flow  case.  Comparisons  of 
the  predicted  and  experimental  forces  and  moment  are  shown  in  Fig.  5.7.  The  standard 
Spalart-Allmaras  model  seems  to  give  the  best  agreement  with  the  experimental  data, 
but  all  the  models  yield  reasonable  results  as  well. 

The  pressure  contours  for  the  airfoil  oscillating  about  the  11°  mean  angle-of-attack 
is  shown  in  Fig.  5.8.  The  airfoil  is  at  a  =  14.93°  on  the  up  stroke.  The  pres¬ 
ence  of  small-scale  structure  in  the  near  wake  of  the  airfoil  can  be  seen  in  the  hybrid 
RANS/LES  results.  Eddy  viscosity  contours  for  the  four  turbulence  models  are  shown 
in  Fig.  5.9.  The  two  RANS  models  predict  similar  distributions  on  the  airfoil,  but  the 
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SST  model  predicts  higher  levels  of  eddy  viscosity  in  the  wake  than  does  the  Spalart- 
Allmaras  model.  Both  of  the  hybrid  RANS/LES  models  are  seen  to  reduce  the  eddy 
viscosity  in  the  aft  region  of  the  airfoil  and  in  the  near  wake.  The  hybrid  SST  model 
recovers  a  higher  level  of  eddy  viscosity  in  the  wake  than  does  the  DES  model.  Com¬ 
parisons  of  the  predicted  and  experimental  forces  and  moment  are  shown  in  Fig.  5.10. 
The  hybrid  RANS/LES  models  and  the  SST  model  show  distinct  oscillations  on  the 
down  stroke  at  the  high  angles-of-attack.  The  Spalart-Allmaras  model  does  not  predict 
this  effect.  The  data  show  some  small  oscillations  in  this  angle-of-attack  range.  It  is 
suspected  that  with  the  limited  number  of  pressure  ports,  the  low  frequency  filtering, 
and  the  three-dimensional  effects  the  experimental  results  have  filtered  and  damped 
out  most  of  the  oscillations.  The  spectral  result  for  the  normal  force  is  shown  in  Fig. 
5.11.  These  results  were  taken  from  the  last  4096  time  steps  of  the  solutions.  The  hy¬ 
brid  RANS/LES  models  and  the  SST  model  all  show  significant  energy  above  500  Hz, 
which  would  not  be  seen  in  the  experimental  results  because  of  the  filters  employed  in 
the  data  reduction. 

The  pressure  contours  for  the  airfoil  oscillating  about  the  15°  mean  angle-of-attack 
is  shown  in  Fig.  5.12.  The  airfoil  is  at  a  =  18.92°  on  the  down  stroke.  The  hybrid 
RANS/LES  models  and  the  SST  model  all  show  a  large  vortical  structure  above  the 
airfoil  with  smaller  structures  coming  off  the  airfoil  leading  and  trailing  edges.  The 
minimum  pressure  for  the  large  vortex  is  lower  for  the  hybrid  RANS/LES  models  than 
in  the  SST  model.  The  Spalart-Allmaras  model  does  not  show  any  structure  above  the 
airfoil.  Eddy  viscosity  contours  for  the  four  models  are  shown  in  Fig.  5.13.  Again 
there  is  considerable  difference  in  the  results.  The  SST  model  has  high  values  of  eddy 
viscosity  in  the  region  above  the  upper  surface  of  the  airfoil.  The  hybrid  SST  model 
has  much  lower  values  near  the  airfoil  surface,  but  approximates  the  SST  model  in  the 
region  away  from  the  airfoil.  The  DES  model  has  low  eddy  viscosity  values  every¬ 
where,  which  again  may  be  due  to  the  strategy  used  to  switch  between  RANS  and  LES 
modes.  The  Spalart-Allmaras  model  has  high  eddy  viscosity  values  near  the  airfoil 
surface.  Comparisons  of  the  predicted  and  experimental  forces  and  moment  are  shown 
in  Fig.  5.14.  All  of  the  models  show  oscillations  on  the  down  stroke  and  on  the  end  of 
the  upstroke.  The  data  does  not  show  this  effect.  Again,  it  is  suspected  that  the  exper¬ 
imental  results  have  filtered  out  most  of  the  oscillations  in  the  data  reduction  process. 
The  spectral  result  for  the  normal  force  is  shown  in  Fig.  5.15.  These  results  were  taken 
from  the  last  4096  time  steps  of  the  solutions.  The  hybrid  RANS/LES  models  show 
significant  energy  above  200  Hz,  which  would  not  be  seen  in  the  experimental  results 
because  of  the  data  reduction  procedure.  The  RANS  models  both  have  a  spectral  peak 
at  about  200  Hz. 


5.3  WICS  Bay 

Unsteady  Navier-Stokes  calculations  were  also  performed  for  the  WICS  (Weapons 
Internal  Carriage  and  Separation)  L/D  =  4.5  weapons  bay  [31]  for  M  =  0.95  and 
Re  =  2.5  x  10 6 / ft .  The  wind  tunnel  data  were  obtained  in  the  AEDC  four-foot  tran¬ 
sonic  wind  tunnel.  The  weapons  bay  was  18  in.  long,  4  in.  wide,  and  4  in.  deep.  The 
computational  geometry  was  a  flat  plate  that  extended  15  in.  upstream  of  the  bay  to 
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match  the  experiment  and  25  in.  downstream  of  the  bay.  The  sides  of  the  computational 
grid  extended  50  in.  on  either  side  of  the  bay  centerline.  As  in  the  previous  simulations, 
wall  functions  for  all  viscous  walls.  The  wall  spacing  was  chosen  as  0.0075  in.,  which 
corresponds  to  a  y+  of  50  on  the  upstream  plate.  The  grid  system  had  1,132,860  points 
broken  into  eight  grids.  A  constant  time  step  of  7.88  x  1 0  6  seconds  was  used  in  the 
calculation.  The  calculation  was  run  2000  steps  to  remove  the  initial  transients  and  all 
unsteady  results  were  processed  over  the  last  8192  time  steps. 

Mach  number  contours  on  the  WICS  bay  centerline  are  shown  in  Fig.  5.16  for  all 
the  turbulence  models.  All  of  the  models  except  the  DES  models  produce  qualitatively 
similar  results.  The  DES  model  seems  to  be  much  more  dynamic  at  this  moment  of 
observation.  The  eddy  viscosity  contours  are  shown  in  Fig.  5.17.  The  Spalart-Allmaras 
model  predicts  high  values  of  eddy  viscosity  in  side  the  bay.  The  SST  model  has 
regions  of  high  eddy  viscosity  within  the  bay,  but  has  less  eddy  viscosity  overall  than 
does  the  Spalart-Allmaras  model.  The  hybrid  SST  model  has  much  lower  values  of 
eddy  viscosity  in  the  bay  than  do  the  RANS  models,  and  the  DES  model  seems  to  have 
almost  no  eddy  viscosity  in  the  bay,  indicating  that  the  simulation  is  virtually  a  MILES 
approach  in  this  region.  Mach  number  contours  at  X/L  =  0.61  are  shown  in  Fig.  5.18. 
These  contours  indicate  the  three-dimensionality  of  the  flow.  The  DES  model  seems  to 
demonstrate  the  largest  three-dimensional  effects. 

The  time-averaged  pressure  coefficient  and  the  sound  pressure  level  along  the  bay 
centerline  are  shown  in  Fig.  5.19.  All  four  of  the  models  are  in  reasonable  agreement 
with  data.  The  three  outlying  sound  pressure  level  data  points  are  suspected  to  be  bad 
data,  but  are  included  here  for  completeness.  Spectral  results  for  two  locations  in  the 
bay  are  shown  in  Fig.  5.20.  All  of  the  models  predict  the  largest  spectral  peak  at  about 
480  Hz.  The  RANS  models  predict  the  first  three  peaks,  while  the  hybrid  RANS/LES 
models  predict  the  first  four  peaks.  The  DES  model  does  seem  to  produce  some  spu¬ 
rious  peaks  not  seen  in  the  data  or  in  any  of  the  other  models.  The  hybrid  RANS/LES 
models  are  in  much  better  agreement  with  the  data  at  the  higher  frequencies  than  are 
the  RANS  models. 


5.4  Planar  Shear  Layer 

One  area  of  great  importance  for  practical  application  of  hybrid  RANS/LES  models  in  a 
production  environment  is  the  examination  of  their  failure  modes.  When,  for  whatever 
reason,  one  has  a  fine  grid  but  little  or  no  unsteadiness  is  reaching  that  part  of  the 
domain,  it  is  important  that  the  turbulence  model  “fail”  in  a  benign  way.  In  the  current 
work,  it  has  been  observed  that  the  DES  model  tends  to  revert  to  a  laminar  solution 
(away  from  walls),  while  the  hybrid  model  falls  back  to  something  closer  to  a  standard 
RANS  solution.  Illustrative  of  this  are  the  results  of  3-D  simulations  using  the  NXAIR 
RANS  flow  solver  of  the  same  mixing  layer  case  discussed  in  Chapter  3.  Plots  of 
mixing  layer  thickness  are  shown  in  Fig.  5.21.  The  coarse  grid  dimensions  in  this  case 
were  188  x  106  x61  and  for  the  fine  grid  the  dimensions  were  188x  106  x  141.  The 
fine  grid  spanwise  spacing  is  half  that  of  the  coarse  grid  at  the  center  of  the  duct.  Both 
cases  were  run  without  any  sort  of  forcing  or  pseudo-turbulence  at  the  inflow  that  would 
trigger  chaotic  3-D  behavior.  For  the  coarse  grid,  all  but  the  DES  model  are  predicting 
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a  reasonable  RANS  mixing  layer  solution.  The  DES  model  is  underpredicting  the 
growth  somewhat.  The  situation  is  even  more  striking  for  the  fine  grid,  however,  where 
the  DES  model  is  clearly  predicting  a  laminar  solution.  In  contrast,  the  hybrid  SST 
model,  although  it  does  show  some  deviation  from  the  RANS  models,  is  still  predicting 
mixing  layer  growth  along  the  lines  of  a  turbulent  RANS  simulation. 

If  this  trend  holds  up  in  further  testing,  then  it  would  represent  a  clear  advantage 
for  the  hybrid  model  over  the  DES  model.  This  is  especially  true  since,  for  complex 
production  applications,  grid  quality  is  often  impossible  to  maintain  over  the  entire 
domain.  As  long  as  sufficient  resolution  is  present  to  capture  the  unsteady  features  of 
interest  in  the  region  of  interest,  it  is  often  acceptable  for  the  rest  of  the  domain  to  revert 
to,  in  effect,  a  turbulent  RANS  solution.  A  reversion  to  laminar  flow,  however,  would 
not  be  acceptable  for  many  applications. 


5.5  Conclusions 

The  hybrid  RANS/LES  models  provided  reasonable  agreement  with  the  data  in  all 
cases,  and  consistently  predicted  higher  unsteady  energy  levels  than  did  the  traditional 
RANS  models.  The  hybrid  models  also  predicted  lower  minimum  pressures  in  vortex 
cores  than  did  the  RANS  models,  which  indicates  that  they  can  be  used  to  remove  a 
known  restriction  of  RANS  type  models. 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.1:  Instantaneous  Pressure  Contours  for  a  Circular  Cylinder  for  M  =  0.2  and 
Red  =  8.0  x  106 
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Figure  5.2:  Pressure  Distribution  Across  the  Nearest  Shed  Vortex  of  a  Circular  Cylinder 
at  M  =  0.2  and  Red  =  8.0  x  106 
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(a)  SST  Model 


(b)  Hybrid  SST  Model 


(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.3:  Instantaneous  Eddy  Viscosity  Contours  (Nondimensionalized  by  Free 
Stream  Molecular  Viscosity)  for  a  Circular  Cylinder  at  M  =  0.2  and  Re,[  =  8.0  x  106 
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Figure  5.4:  PSD  for  the  Circular  Cylinder  Axial  Force  at  M  =  0.2  and  Red  =  8.0  x  106 


Strouhal  Number 


Figure  5.5:  PSD  for  the  Circular  Cylinder  Normal  Force  at  M  =  0.2  and  Red  =  8.0  x 
106 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.6:  Pressure  Contours  on  a  NACA  0015  airfoil  at  a  =  7.92°  on  the  Down 
Stroke  Oscillating  About  ao  =  4°  with  K  =  0.1 
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(a)  Drag  Coeffi  cient 


(b)  Lift  Coeffi  cient 


(c)  Moment  Coeffi  cient 


Figure  5.7:  Forces  and  Moments  on  a  NACA  0015  Airfoil  Oscillating  About  oto  =  4° 
With  a  Reduced  Frequency  of  K  =  0.1 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.8:  Pressure  Contours  on  a  NACA  0015  airfoil  at  a  =  14.93°  on  the  Up  Stroke 
Oscillating  About  oto  =  11°  with  K  =  0.1 
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(a)  SST  Model 


(b)  Hybrid  SST  Model 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.9:  Eddy  Viscosity  Contours  on  a  NACA0015  airfoil  at  a  =  14.93°  on  the  Up 
Stroke  Oscillating  About  ao  =  11°  with  K  =  0. 1 
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(a)  Drag  Coeffi  cient 


(b)  Lift  Coeffi  cient 


(c)  Moment  Coeffi  cient 


Figure  5.10:  Forces  and  Moments  on  a  NACA  0015  Airfoil  Oscillating  About  ao  =  1 1  ° 
With  a  Reduced  Frequency  of  K  =  0.1 
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Figure  5.11:  Predicted  PSD  for  the  Normal  Force  on  a  NACA  0015  Airfoil  Oscillating 
About  ao  =  1 1  °  With  a  Reduced  Frequency  of  K  =  0. 1 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.12:  Pressure  Contours  on  a  NACA  0015  airfoil  at  a  =  18.92°  on  the  Down 
Stroke  Oscillating  About  ao  =  15°  with  K  =  0. 1 
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(a)  SST  Model 


(b)  Hybrid  SST  Model 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.13:  Eddy  Viscosity  Contours  on  a  NACA  0015  airfoil  at  a  =  18.92°  on  the 
Down  Stroke  Oscillating  About  ao  =  15°  with  K  =  0. 1 
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(a)  Drag  Coffi  cient 


(b)  Lift  Coeffi  cient 


(c)  Moment  Coeffi  cient 


Figure  5.14:  Forces  and  Moments  on  aNACA  0015  Airfoil  Oscillating  About  ao  =  15° 
With  a  Reduced  Frequency  of  K  =  0.1 
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Figure  5.15:  Predicted  PSD  for  the  Normal  Force  on  a  NACA  0015  Airfoil  Oscillating 
About  ao  =  15°  With  a  Reduced  Frequency  of  K  =  0. 1 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.16:  Instantaneous  Mach  Number  Contours  on  the  WICS  Bay  Centerline  for 
M  =  0.95 
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(b)  Hybrid  SST  Model 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.17:  Instantaneous  Eddy  Viscosity  Contours  on  the  WICS  Bay  Centerline  for 
M  =  0.95 
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(c)  Spalart-Allmaras  Model 


(d)  DES  Model 


Figure  5.18:  Instantaneous  Mach  Number  Contours  at  x/D  =  0.6  for  the  WICS  Bay  at 
M  =  0.95 
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a)  Front  Wall 
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c)  Back  Wall 


Figure  5.19:  Time  Averaged  Pressure  Coefficient  and  Sound  Pressure  Level  on  the 
WICS  Bay  Centerline 
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(a)  WICS  bay  fbor  centerline  X/L  =  0.985 


(b)  WICS  bay  back  wall  centerline  at  Z/H  =  0.82 


Figure  5.20:  FFT  for  two  locations  in  the  WICS  bay  at  M  =  0.95 
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SHEAR  LAYER  COARSE  GRID  Ml  =0.51  M2=1.8 


(a)  Coarse  Grid 


SHEAR  LAYER  FINE  GRID  M1=0.51  M2=1.8 


(b)  Fine  Grid 


Figure  5.21:  Mixing  Layer  Growth  Rate  Predicted  Using  a  RANS  Flow  Solver  With 
No  Inflow  Perturbations 
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Chapter  6 

Code  Veri  cation  Using  the 
Method  of  Manufactured 
Solutions 

6.1  Introduction 

As  mentioned  earlier,  technology  transition  has  been  a  major  component  of  this  project 
from  the  beginning.  It  is  not  sufficient  to  merely  run  simplified  cases  using  a  research 
code;  The  resulting  technology  must  be  transitioned  to  a  production  flow  solver  if  at  all 
possible.  As  discussed  above,  the  flow  solver  chosen  for  this  effort  was  the  Wind-US 
code.  Wind-US  is  a  relatively  large  and  complex  code  with  numerous  options  for  al¬ 
most  every  aspect  of  a  CFD  simulation.  As  is  characteristic  of  complex  systems,  it  has 
in  the  past  been  found  that  updates  to  fix  one  problem  often  cause  other  problems.  It 
therefore  becomes  very  important  to  verify  that  the  algorithms  required  to  perform  the 
Facility  Unsteadiness  Research  project  remain  functional  through  all  the  code  modi¬ 
fications.  Indeed,  much  time  was  spent  in  FY2004  tracking  down  a  problem  which 
resulted  in  Newton  iterations  not  working  properly  in  combination  with  the  second  or¬ 
der  implicit  time  marching  scheme.  The  Method  of  Manufactured  Solutions  has  been 
found  to  be  extremely  useful  in  this  regard,  and  therefore  this  section  covering  it  is 
included  in  this  final  report. 

6.2  Overview  of  the  Method  of  Manufactured  Solutions 

The  analysis  of  error  in  computational  fluid  dynamics  (CFD)  codes  can  be  divided 
into  two  general  categories;  verification  and  validation.  Validation  refers  to  whether 
the  correct  equations  are  solved  for  the  problem  at  hand.  For  example,  a  linear  Euler 
solver  will  not  be  able  to  accurately  resolve  strong  shocks,  no  matter  how  well  im¬ 
plemented.  Verification  addresses  the  issue  of  whether  or  not  the  chosen  governing 
equations  are  solved  correctly.  In  verifying  CFD  results,  one  must  address  such  things 
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as  the  sufficiency  of  grid  resolution,  completeness  of  convergence  (or  accuracy  of  time 
marching),  and  the  presence  (or  not)  of  coding  errors. 

One  tool  for  code  verification  which  has  only  recently  begun  to  be  applied  to  large, 
general-purpose  flow  solvers  is  the  method  of  manufactured  solutions  (MMS).  This 
technique  is  a  rigorous  method  for  finding  and  eliminating  coding  mistakes.  While 
any  given  series  of  MMS  runs  can  only  examine  a  subset  of  the  available  options  in  a 
modern  production  flow  solver,  within  that  subset,  the  results  that  it  yields  are  fairly 
conclusive.  MMS  was  originally  developed  by  Roache  and  Steinburg[32]  and  extended 
by  Roache  et  al.  [33] .  The  term  “manufactured”  solutions  was  coined  by  Oberkampf 
and  Blottner[34]  and  refers  to  the  fact  that  the  solutions  are  arbitrarily  chosen  (or  man¬ 
ufactured).  The  governing  equations  are  then  modified  by  the  addition  of  source  terms 
such  that  the  manufactured  solution  satisfies  the  governing  equations.  The  work  of 
Salari  and  Knupp[35]  contains  a  detailed  discussion  of  MMS  as  applied  to  a  variety  of 
partial  differential  equations. 

The  actual  application  of  the  method  of  manufactured  solutions  is  now  discussed 
in  somewhat  more  detail.  The  overall  procedure  is  as  follows: 

•  Choose  the  form  of  the  governing  equations 

•  Chose  the  form  of  the  manufactured  solution 

•  Apply  the  governing  equations  to  the  manufactured  solution  to  generate  analyti¬ 
cal  source  terms 

•  Solve  the  equations  on  multiple  mesh  levels  using  the  source  terms 

•  Evaluate  the  global  discretization  error  in  the  numerical  solutions 

•  Determine  the  order  of  accuracy 

The  analytical  solutions  can  be  almost  anything,  but  to  be  useful,  they  must  at  least 
satisfy  the  following  conditions: 

•  They  must  be  continuous.  Current  numerical  schemes  all  reduce  to  first  order  at 
discontinuities-if  they  can  handle  them  at  all. 

•  Each  variable  in  the  solution  must  be  continuously  differentiable  up  to  the  order 
of  the  corresponding  terms  in  the  governing  equations.  One  can  only  test  the 
terms  of  the  governing  equations  for  which  the  solutions  give  non-zero  contribu¬ 
tions 

•  The  Taylor  series  expansion  of  the  solutions  must  include  non-zero  terms  at  least 
up  to  the  purported  order  of  the  numerical  scheme  to  be  tested.  As  above,  only 
non-zero  contributions  from  the  manufactured  solutions  can  be  tested. 

•  In  addition,  for  best  effect,  the  contribution  from  each  term  in  the  governing 
equations  should  be  of  the  same  order  of  magnitude.  This  prevents  the  larger 
terms  from  masking  errors  in  other  terms  of  smaller  magnitude.  For  example, 
in  the  current  work,  when  the  Navier-Stokes  equations  are  solved,  the  viscosity 
is  fixed  at  a  value  of  10  A  ■  s/m2 in  order  to  bring  the  contribution  of  the  viscous 
terms  up  to  the  level  of  the  inviscid  terms. 
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•  In  order  to  avoid  numerical  problems,  the  solutions  should  avoid  negative  values 
of  quantities  which  should  physically  only  be  positive  (e.g.  pressure,  density, 
and  temperature). 

For  the  cases  discussed  here,  the  analytical  solutions  were  chosen  to  be  functions  of 
sines  and  cosines  designed  such  that  each  term  in  the  governing  equations  to  be  exam¬ 
ined  (either  the  Euler  equations  or  the  laminar  Navier-Stokes  equations)  would  produce 
a  non-zero  value.  For  the  Euler  equations,  the  manufactured  solution  for  each  specified 
variable  was  of  the  form: 


§Euler  (v,  Y. 7. )  —  <j)0  +  tyxfyx 


(a$xx  +  cx)  K 
L 

/ 


-$yf b 


,yy+cy)  ns 


(6.1) 


where  <|)  =  p,  u,  v,  w,  or  p  and  the  /  functions  denote  either  the  sine  or  cosine  func¬ 
tion.  L  is  a  reference  length  (taken  as  1  m  here).  The  various  <f>’s,  c/’s,  and  c’s  that 
appear  on  the  right  hand  side  are  constant  coefficients.  Note  that  subscripts  do  not 
denote  differentiation  in  this  case.  The  constants  used  in  the  various  tests  of  the  Euler 
equations  solver  are  contained  in  Tables  6.1,  6.2,  6.3,  and  6.4. 

For  the  Navier-Stokes  solutions,  the  solutions  take  the  form  of  the  Euler  equations 
solutions  with  additional  cross  terms  added  to  ensure  that  the  viscous  terms  are  fully 
exercised: 

§NS  (x,y,z)  =  <\>Euler(x,y,z)+  tyxyfyty 


The  constants  used  for  the  Navier-Stokes  equations  tests  are  shown  in  Tables  6.5,  6.6, 
6.7,  and  6.8. 

A  symbolic  manipulation  software  package  was  used  to  analytically  differentiate 
the  general  manufactured  solution  (either  equation  6.1  or  6.2)  according  to  the  govern¬ 
ing  equations.  Since  the  manufactured  solutions  do  not  actually  solve  the  governing 
equations,  this  will  result  in  extra  terms  which  are  then  added  as  source  terms  to  drive 
the  solution  toward  the  specified  manufactured  solution.  As  the  source  terms  are  quite 
complex,  space  considerations  prevent  them  being  included  in  this  report. 


6.3  Numerical  Scheme 

Wind-US  contains,  among  other  things,  a  structured  compressible  Navier-Stokes  solver. 
It  uses  a  finite  volume  discretization  in  the  solution  of  the  governing  equations.  Inte¬ 
rior  mesh  points  are  treated  as  cell  centers  and  the  corners  that  define  a  given  control 
volume  around  them  are  computed  by  averaging  the  coordinates  of  neighboring  cell 
centers.  With  this  information,  cell  volumes  and  face  areas  can  be  readily  computed. 
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On  the  boundaries,  however,  Wind-US  uses  a  technique  that  was  inherited  from  the 
NASTD  code.  Rather  than  solving  half  cells  at  boundaries,  the  nearest  interior  cell  is 
enlarged  to  reach  all  the  way  to  the  boundary.  Therefore,  if  a  uniform  grid  is  input,  this 
results  in  so-called  “fat”  cells  at  boundaries,  as  illustrated  in  Fig  6.1.  These  fat  cells 
have  two  main  drawbacks:  first,  they  introduce  a  discontinuity  in  in  the  grid  spacing  at 
the  boundaries,  and  second,  the  nominal  “center”  of  a  fat  cell  is  not  actually  at  the  cen¬ 
ter  of  mass  of  the  volume  (even  on  a  uniform  grid).  Both  of  these  issues  are  potential 
sources  of  error.  In  order  to  avoid  these  problems  in  the  current  work,  a  preproces¬ 
sor  has  been  employed  which  reads  in  a  smooth  well-defined  nodal  grid,  computes 
cell  centers,  and  then  adds  a  fringe  of  boundary  points  around  the  edges.  When  such 
a  grid  is  loaded  into  Wind-US,  its  internal  scheme  will  create  a  more  uniform  set  of 
finite  volumes  (see  Fig.  6.2).  For  a  truly  uniform  Cartesian  nodal  grid,  the  resulting 
“cell-centered"  grid  has  uniform  volume  and  the  cell  centers  are  exactly  in  the  center  of 
mass.  Thus,  Wind-US  behaves  much  like  a  cell-centered  finite-volume  code  although 
the  internal  algorithm  remains  node-centered. 

By  default  Wind-US  uses  a  Roe  scheme  on  the  right  hand  side  with  a  flux  extrap¬ 
olation  algorithm  that  attempts  to  account  for  grid  stretching  to  maintain  its  order  on 
non-uniform  meshes.  For  the  Navier-Stokes  results,  the  viscous  fluxes  are  evaluated  at 
cell  faces  using  a  central-difference  scheme.  This  results  in  an  overall  scheme  which, 
on  uniform  meshes  is  formally  second  order  accurate  in  space.  While  it  can  be  run 
explicitly,  Wind-US  is  by  default  an  implicit  code,  and  the  standard  inversion  method 
is  an  approximate  factorization  scheme.  Although  there  are  numerous  alternatives  for 
both  the  right  and  left  hand  side,  the  discussion  here  will  attempt  to  focus  on  the  code 
options  of  interest  when  running  large  scale  unsteady  simulations.  Additional  informa¬ 
tion  about  the  numerics  of  the  structured  portion  of  the  Wind-US  code  may  be  found 
in  the  work  of  Bush[36],  Bush  et  al. [37],  and  Power  and  Underwood[38]. 

6.4  Test  Cases 

MMS  has  been  used  to  test  the  Wind-US  flow  solver  in  several  configurations.  The 
first  set  of  cases  were  run  on  a  two-dimensional  uniform  Cartesian  mesh  using  the 
default  numerical  scheme.  For  the  two-dimensional  cases,  results  were  obtained  (or 
attempted)  for  seven  resolutions  ranging  from  8  x  8  to  128  x  128.  The  domain  for  all 
these  cases  was  a  one  meter  square.  Previous  applications  of  MMS  to  the  Wind  solver 
for  2-D  supersonic  Euler  and  subsonic  Navier-Stokes  flows  have  been  performed  by 
the  Roy  et  al.[39].  For  the  current  work,  multi-block  versions  of  the  same  cases  were 
run  to  assess  the  effect  of  block-to-block  coupling  on  the  scheme’s  overall  order  of 
accuracy.  Also,  for  this  work,  subsonic  Euler  and  supersonic  Navier-Stokes  cases  are 
investigated.  Similar  runs  are  also  performed  on  3-D  versions  of  the  above  cases.  For 
the  3-D  runs,  however,  machine  limitations  prevented  the  running  of  96  x  96  x  96  and 
128  X  128  x  128  single  block  grids.  A  multi-block  version  of  the  96  x  96  x  96  case  was 
generated  using  eight  48  x  48  x  48  grids. 

Of  particular  interest  for  this  project  are  the  higher  order  algorithms  available  in 
the  code.  Thus,  the  third,  fourth,  and  fifth  order  schemes  for  computing  inviscid  fluxes 
have  been  tested. 
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Finally,  the  order  of  accuracy  of  the  Wind-US  code  on  curvilinear  grids  is  investi¬ 
gated.  The  grid  used  for  the  curvilinear  test  case  was  originally  designed  for  an  internal 
duct  flow.  It  features  clustering  and  curvature  to  some  degree  in  all  three  coordinate 
directions. 

For  all  of  the  subsonic  cases,  outer  boundaries  were  specified  to  match  the  pre¬ 
scribed  manufactured  solutions.  For  the  supersonic  cases,  however,  outflow  boundaries 
were  solved  using  the  code’s  normal  supersonic  outflow  boundary  condition  (which 
simply  extrapolates  variables  from  the  interior)  unless  otherwise  noted. 


6.5  Results 

Results  for  the  2-D  supersonic  Euler  case  using  the  default  numerical  scheme  of  Wind- 
US  are  shown  in  Fig.  6.3.  The  Lt norm  of  the  error  at  convergence  is  plotted  for  the 
conserved  variables  at  various  grid  resolutions.  The  results  for  single  block  grids  are 
shown  with  open  symbols,  while  the  filled  symbols  represent  a  multi-block  run  (the 
96  x  96  grid  was  split  into  four  equal  parts).  At  the  higher  resolutions  (48  x  48  and 
above),  the  code  is  clearly  showing  second  order  convergence.  The  “bulge”  in  the 
curve  for  the  32  x  32  grid  is  caused  by  the  flux  limiter  activating  at  this  level,  which 
reduces  the  scheme  to  first  order.  The  proof  of  this  is  found  in  Fig.  6.4,  which  has 
results  from  the  same  case,  but  this  time  the  flux  limiter  was  deactivated.  For  this 
case,  the  code  shows  second  order  convergence  throughout  the  range  of  grids,  although 
for  the  coarsest  grid  (with  8x8  resolution)  the  convergence  seems  to  be  dropping  off 
somewhat.  Notable  on  both  of  these  figures,  the  converged  error  of  the  multi-block  case 
is  almost  identical  to  the  corresponding  single  block  case.  The  differences  that  exist 
are  due  to  the  implicit  scheme  not  having  information  available  from  the  entire  domain 
when  it  solves  each  block  in  the  multi-block  case.  The  equivalent  three-dimensional 
case  produced,  as  expected,  much  the  same  behavior  as  the  two-dimensional  results 
(see  Fig.  6.5). 

Results  from  the  subsonic  Euler  case  are  shown  in  Fig.  6.6  for  two  dimensions  and 
Fig.  6.7  for  three.  For  both  these  cases,  the  code  clearly  demonstrates  second  order 
convergence.  Again,  the  equivalent  multi-block  solutions  also  appear  to  be  behaving 
much  like  their  single-block  counterparts. 

The  first  Navier-Stokes  case  to  be  examined  is  the  subsonic  case.  Figure  6.8  shows 
the  results  from  the  2-D  case.  The  convergence  rate  at  the  coarser  resolutions  appears  to 
be  somewhat  less  than  second  order  (but  better  than  first  order).  At  higher  resolutions, 
however,  the  curve  approaches  a  second  order  slope,  indicating  that  with  sufficient 
resolution,  the  theoretical  order  can  be  achieved.  The  results  from  the  equivalent  3- 
D  case  are  shown  in  Fig.  6.9.  As  with  the  2-D  case,  the  code  shows  clear  second 
order  behavior.  For  both  of  these  cases,  the  multi-block  runs  failed  to  complete.  This  is 
because  the  algorithm  that  Wind-US  uses  to  communicate  between  abutting  boundaries 
currently  only  passes  the  inviscid  fluxes.  When  the  viscous  fluxes  are  of  the  same 
magnitude  as  the  inviscid  fluxes  (as  in  this  case),  the  coupling  mechanism  breaks  down. 
It  should  be  noted  that,  in  theory,  the  overlapped  coupling  in  Wind-US  has  no  such 
limitations,  but  this  has  yet  to  be  tested  using  MMS. 

The  power  of  MMS  as  a  verification  tool  was  highlighted  by  the  results  from  the 
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2-D  supersonic  Navier-Stokes  test  case.  Figure  6.10  the  convergence  of  the  origi¬ 
nal,  default  scheme  with  first  order  extrapolation  at  outflow  boundaries.  As  the  fig¬ 
ure  makes  painfully  clear,  Wind-US  was  not  behaving  as  a  second  order  code  at  all. 
Since  the  Euler  cases  already  verified  the  inviscid  fluxes  as  being  correct,  and  the  sub¬ 
sonic  Navier-Stokes  case  seemed  to  indicate  that  the  viscous  terms  were  also  correct, 
suspicion  quickly  focused  on  the  outflow  boundary  condition  (recall  that  the  subsonic 
case  specified  the  exact  solution  on  outflow  boundaries  where  the  supersonic  case  uses 
Wind-US’s  normal  outflow  boundary  condition).  To  test  this,  the  case  was  re-run  using 
zeroth  order  extrapolation  on  the  outflow. 

As  seen  on  Fig.  6.11,  these  results  are  markedly  different  than  the  previous,  and 
indeed  show  something  of  a  zeroth  order  rate  of  convergence  at  the  higher  resolutions. 
This  indicates  that  not  only  was  the  outflow  boundary  affecting  the  solution  (via  the 
viscous  terms),  but  that  at  higher  resolutions  it  was  the  dominant  factor.  Upon  exam¬ 
ination,  the  extrapolation  algorithm  was  found  to  have  two  significant  bugs.  First  of 
all,  the  basic  algorithm  was  inappropriate  for  the  pseudo-cell-centered  grid  which  was 
being  employed.  Second,  an  attempt  was  being  made  to  apply  a  flux  limiter  to  the 
extrapolation,  and  this  algorithm  was  completely  unsuitable  for  the  task.  Having  fixed 
these  two  problems,  the  case  was  re-run  with  first  order  extrapolation  on  the  bound¬ 
aries.  Figure  6.12  shows  that  these  changes  resulted  in  a  dramatic  improvement  in  the 
error.  The  code  now  clearly  shows  roughly  second  order  convergence  on  the  coarse 
grids,  while  the  fine  grids  have  a  first  order  convergence  rate. 

Seeing  that  the  Euler  case  had  been  affected  by  the  presence  of  the  flux  limiter,  it 
was  decided  to  re-run  these  cases  without  it.  As  Fig.  6.13  shows,  there  was  little  or  no 
difference  in  the  results.  This  led  to  the  conclusion  that  the  outflow  boundaries  must  be 
the  source  of  the  first  order  behavior.  To  prove  that  the  influence  of  the  boundary  was 
driving  the  first  order  results,  a  final  set  of  runs  was  made  which  specified  the  analytic 
solution  on  all  boundaries.  The  results  from  this  case  (shown  in  Fig.  6.14)  strongly  in¬ 
dicate  that  it  is  in  fact  the  first  order  nature  of  the  boundary  condition  which  was  driving 
the  previously  observed  behavior.  The  slight  glitch  at  the  highest  resolutions  resisted 
explanation  until  very  recently,  when  this  case  was  rerun  using  the  Koren  limiter.  With 
the  Koren  limiter  engaged,  the  results  (in  Fig.  6.15)  show  fully  second  order  behavior 
throughout  the  range  of  grid.  It  appears  that  the  bump  that  appears  when  either  the 
default  TVD  limiter  or  no  limiter  at  all  is  used  is  due  to  error  building  up  in  the  domain 
that  the  scheme  cannot  naturally  dissipate.  The  default  Superbee  limiter,  which  should 
counter  this  noise  is,  it  appears,  capable  of  introducing  noise  of  its  own  into  the  solu¬ 
tion  (counter  to  its  design).  An  equivalent  set  of  three-dimensional  grid  runs  was  made 
using  the  first  order  boundary  condition.  As  one  would  expect,  the  results  (shown  in 
Fig.  6.16)  are  very  much  in  line  with  the  2-D  case. 

The  multi-block  results  for  the  supersonic  Navier-Stokes  cases  were  somewhat 
mixed.  With  the  buggy  first  order  or  the  zeroth  order  boundary  conditions,  the  multi¬ 
block  results  were  not  too  different  from  the  single  block.  As  the  boundary  condition’s 
effects  are  reduced,  however,  the  problems  with  the  inviscid-only  block  coupling  algo¬ 
rithm  become  more  evident.  Interestingly,  the  3-D  run  did  not  show  as  much  variation, 
although  with  only  one  data  point,  nothing  is  conclusive.  Up  until  now,  only  abut¬ 
ting,  point-matched  boundaries  have  been  considered.  Wind-US  is  also  capable  of 
computing  overlapped  boundary  conditions  with  fringes  of  arbitrary  size.  In  theory. 
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a  point-matched  double  fringe  overlapped  zonal  coupling  configuration  ought  to  en¬ 
able  fully  second  order  behavior  throughout  the  domain.  To  test  this,  the  3-D  subsonic 
Navier-Stokes  case  was  re-run  after  splitting  the  grids  into  eight  blocks  each.  As  the 
results  in  Fig.  6.17  show,  second  order  behavior  was  predicted  in  this  configuration. 
Note  that  for  this  case,  it  was  necessary  to  run  with  “fat  cells”  at  the  boundaries  in 
order  to  ensure  that  the  overlapping  regions  were  point-matched.  This  results  in  an 
increased  error,  as  can  be  seen  by  comparing  with  Fig.  6.9.  The  Lin orm  of  the  error  is, 
on  average,  increased  by  roughly  a  factor  of  two  for  the  case  with  “fat  cells”. 

For  the  facility  unsteadiness  work,  it  would  be  desirable  to  have  the  option  of  run¬ 
ning  even  higher  order  schemes  (at  least  in  regions  where  the  flow  is  smooth).  There¬ 
fore,  the  third,  fourth,  and  fifth  order  upwind-biased  schemes  were  investigated.  The 
first  results  obtained  with  the  third  order  scheme  (in  Fig.  6.18)  show,  surprisingly,  a 
second-order  convergence  pattern.  Suspecting  more  limiter  problems,  the  case  was 
rerun  without  a  TVD  limiter.  The  results,  shown  in  Fig.  6.19,  indicate  that  the  base 
scheme  is  correctly  third  order.  The  question  arose,  however,  as  to  whether  or  not 
there  was  something  inherently  second  order  about  the  way  limiters  are  implemented 
in  Wind-US.  To  test  this,  the  Koren  limiter  was  used  for  another  sequence  of  runs.  This 
case  showed  (Fig.  6.20)  that  the  problems  were  due  to  something  specific  to  the  default 
limiter,  and  that,  provided  the  solution  was  smooth,  the  Koren  limiter  allows  third  order 
behavior. 

For  the  fourth  and  fifth  order  schemes,  no  TVD  limiting  is  available,  so  none  of  the 
problems  seen  with  the  second  and  third  order  schemes  were  expected.  As  shown  in 
Fig.  6.21,  however,  while  it  appears  that  the  coarsest  grids  show  approximate  fourth 
order  convergence,  a  strongly  second  order  behavior  is  seen  on  the  finer  grids.  For  the 
fifth  order  scheme,  the  second  order  behavior  is  even  more  dominant  (see  Fig.  6.22). 
The  preliminary  analysis  of  this  problem  indicated  that  it  was  due  to  a  lack  of  numerical 
precision.  The  error  level  for  these  schemes  is  a  factor  of  two  smaller  than  for  the  third 
order  scheme,  and  an  order  of  magnitude  less  than  the  default  second  order  scheme. 
Also,  the  behavior  at  the  finer  resolutions  is  very  similar  for  both  the  fourth  and  fifth 
order  schemes.  Some  efforts  were  made  to  increase  the  variability  of  the  manufactured 
solution,  in  order  to  increase  the  error  magnitude,  but  only  the  coarsest  grids  were 
affected.  This  indicated  that  numerical  precision  was  at  issue,  but  this  could  not  be 
fully  tested  until  a  double  precision  version  of  the  code  became  available  FY2004. 

When  the  fourth  and  fifth  order  cases  were  re-run  using  a  double  precision  version 
of  Wind-US,  it  was  found  that  the  solution  error  did  indeed  show  some  effects  from 
precision.  With  the  double  precision  code,  it  was  possible  to  converge  the  solution  at 
least  a  further  six  orders  of  magnitude.  The  enabled  a  clear  separation  of  the  residual 
error  (due  to  lack  of  convergence)  and  the  error  due  to  the  limitations  of  the  numerical 
scheme.  Unfortunately,  the  overall  error  levels  remained  virtually  unchanged.  The 
source  of  the  second  order  behavior  in  the  higher  order  schemes,  therefore,  remains  a 
mystery. 

Finally,  some  preliminary  MMS  runs  have  been  made  with  a  curvilinear  grid.  The 
grid  (and  converged  density  contours)  are  shown  in  Fig.  6.23.  This  grid  was  originally 
used  for  duct  simulations.  The  topology  is  that  of  an  H-C  grid,  with  a  singular  axis 
running  down  the  middle  of  the  duct.  As  Fig.  6.24  shows,  on  this  grid,  the  convergence 
appears  to  fall  somewhere  between  first  order  and  second.  This  is  probably  the  best 
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that  can  be  hoped  for,  since  the  physical  space-based  extrapolation  does  not  take  grid 
curvature  into  account. 


6.6  Conclusions 

The  method  of  manufactured  solutions  has  thus  far  been  employed  to  verify  the  theo¬ 
retical  order  of  accuracy  for  several  of  the  algorithms  available  in  the  Wind-US  flow 
solver.  Based  on  the  results  presented  above,  one  can  conclude  that,  absent  intervention 
from  the  flux  limiter  or  the  build-up  of  spurious  noise  in  the  solution,  the  default  inte¬ 
rior  scheme  of  Wind-US  has  been  shown  to  be  second  order  accurate  on  uniform  grids 
for  inviscid  flow.  Likewise,  the  third  order  scheme  is  clearly  correctly  implemented, 
although  in  that  case,  the  Koren  limiter  must  be  used  if  TVD  limiting  is  required.  For 
the  fourth  and  fifth  order  schemes,  it  appears  that  a  double  precision  version  of  the 
code  will  be  required  to  verify  the  order  of  accuracy.  There  are  also  strong  indications 
that,  for  inviscid  flow,  the  block-to-block  coupling  scheme  is  also  second  order  (as 
advertised),  at  least  for  point-matched  grids. 

For  the  viscous  cases,  it  has  been  shown  that  the  interior  viscous  scheme  is  second 
order.  The  block  coupling  algorithm,  as  expected,  does  not  perform  well  for  abutting 
grids  when  significant  viscous  terms  are  present. 

The  experience  with  the  supersonic  Navier-Stokes  case  clearly  shows  the  power  of 
MMS  to  identify  coding  problems  and  verify  that  they  have  been  fixed.  It  also  served 
as  a  reminder  that,  no  matter  what  the  nominal  order  of  the  scheme,  the  results  cannot 
exceed  the  accuracy  of  the  boundary  conditions  when  those  boundary  conditions  are 
capable  of  influencing  the  interior  grid. 

Results  from  general  curvilinear  grids  indicate  that,  on  a  smooth  grid  with  reason¬ 
able  grid  resolution,  the  default  scheme  while  not  fully  second  order  accurate,  is  at  least 
better  than  first  order. 
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Variable 

<1)0 

4 

Pa 

«<t>* 

Cx 

f K 

Pv 

a§y 

cy 

P  {kg/m3) 

1 

sin 

0.15 

1 

1 

3 

COS 

-0.1 

0.5 

0.5 

u  {m/s) 

70 

sin 

5 

1.5 

l 

cos 

-7 

0.6 

0 

t'  {m/s) 

90 

cos 

-15 

0.5 

1.5 

sin 

8.5 

2 

3 

1 

p  (N /m1) 

1  x  10s 

cos 

0.2  x  10s 

2 

0 

sin 

0.5  x  10s 

i 

1.4 

Table  6.1:  Manufactured  solution  coefficients  for  2-D  subsonic  Euler  cases 


Variable 

Po 

4 

Pa- 

«<!>, 

Cx 

4 

Pv 

«< |>v 

cy 

w  {kg/m3) 

80 

sin 

10 

1 

T 

0 

sin 

2 

1.5 

1.25 

Variable 

4 

cz 

P  {kg/m3) 

sin 

-0.12 

1.5 

0 

u  (m/s) 

cos 

-1.8 

0.5 

0.5 

v  (m/s) 

sin 

-3 

1.25 

0.75 

w  (m/s) 

cos 

3.5 

1 

0.25 

p  ( N/m -) 

cos 

-0.35  x  10s 

1 

3 

0.75 

Table  6.2:  Additional  Manufactured  solution  coefficients  for  3-D  subsonic  Euler  cases 


Variable 

po 

4 

Pa 

«4>a 

Cx 

4 

Pv 

a$)y 

Cy 

P  (kg/m3) 

1 

sin 

0.15 

1 

1 

3 

COS 

-0.1 

0.5 

0.5 

u  (m/s) 

800 

sin 

50 

1.5 

l 

cos 

-30 

0.6 

0 

t'  (m/s) 

800 

cos 

-75 

0.5 

1.5 

sin 

40 

2 

3 

1 

p  (N/m1) 

1  x  10s 

cos 

0.2  x  10s 

2 

0 

sin 

0.5  x  10s 

i 

1.4 

Table  6.3:  Manufactured  solution  coefficients  for  2-D  supersonic  Euler  cases 


Variable 

Po 

4 

Pa 

cx 

4 

Pv 

Cy 

w  (kg/m3) 

800 

sin 

15 

1 

I 

0 

sin 

-25 

1.5 

1.25 

Variable 

4 

C<\>z 

cz 

P  (kg/m3) 

sin 

-0.12 

1.5 

0 

u  (m/s) 

cos 

-18 

0.5 

0.5 

v  (m/s) 

sin 

-30 

1.25 

0.75 

w  (m/s) 

cos 

35 

1 

0.25 

p  (N/m1) 

cos 

-0.35  x  10s 

1 

T 

0.75 

Table  6.4:  Additional  Manufactured  solution  coefficients  for  3-D  supersonic  Euler 
cases 
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Variable 

<l>0 

4 

4>* 

Cx 

4 

4>y 

a<Dv 

cy 

P  (kg/m3) 

1 

sin 

0.1 

0.75 

0 

cos 

0.15 

1 

0 

u  (m/s) 

70 

sin 

4 

5 

3 

0 

cos 

-12 

1.5 

0 

v  (m/s) 

90 

cos 

-20 

1.5 

0 

sin 

4 

1 

0 

p  ( N/m 2) 

1  x  105 

cos 

-0.3  x  105 

1 

0 

sin 

0.2  x  105 

1.25 

0 

Variable 

K 

<1 by 

P  {kg/m*) 

COS 

0.08 

1.25 

u  (m/s) 

cos 

7 

0.6 

v  (m/s) 

cos 

-11 

0.9 

p  (N/m1) 

sin 

-0.25  x  105 

0.75 

Table  6.5:  Manufactured  solution  coefficients  for2-D  subsonic  Navier-Stokes  cases 


Variable 

<|>o 

A, 

cx 

Ay 

a§y 

Cy 

Uxy 

a§xy 

w  ( kg/ nr !) 

80 

sin 

10 

l 

3 

0 

sin 

2 

1.5 

0 

sin 

8 

0.4 

Variable 

/<t>z 

a$z 

cz 

ftyxz 

<hz 

Uyz 

fyyz 

ayz 

P  (kg/m3) 

sin 

-0.12 

1.5 

0 

sin 

0.09 

1.5 

COS 

0.09 

1 

u  (m/s) 

cos 

-1.8 

0.5 

0 

sin 

6 

0.8 

sin 

6 

0.4 

v  (m/s) 

sin 

-3 

1.25 

0 

cos 

12 

0.7 

sin 

12 

0.8 

w  (m/s) 

cos 

3.5 

1 

0 

sin 

-8 

0.6 

cos 

-8 

0.8 

p  (V/m2) 

cos 

-0.35  x  105 

l 

3 

0 

cos 

0.3  x  105 

1 

cos 

0.3  x  105 

0.5 

Table  6.6:  Additional  Manufactured  solution  coefficients  for  3-D  subsonic  Navier- 
Stokes  cases 


Variable 

<l>0 

4 

4>* 

@§x 

Cx 

4 

4>y 

cy 

P  {kg/m3) 

1 

sin 

0.1 

0.75 

0 

COS 

0.15 

1 

0 

u  (m/s) 

800 

sin 

50 

5 

3 

0 

cos 

-30 

1.5 

0 

v  (m/s) 

800 

cos 

-75 

1.5 

0 

sin 

40 

1 

0 

p  (N/m1) 

1  x  105 

cos 

-0.3  x  105 

1 

0 

sin 

0.2  x  105 

1.25 

0 

Variable 

4. 

§xy 

a§xy 

P  (kg/m3) 

COS 

0.08 

1.25 

u  (m/s) 

cos 

70 

0.6 

v  (m/s) 

cos 

-110 

0.9 

p  (N/m1) 

sin 

-0.25  x  105 

0.75 

Table  6.7:  Manufactured  solution  coefficients  for  2-D  supersonic  Navier-Stokes  cases 
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Variable 

4>o 

4 

<j>* 

a§x 

Cx 

4 

ah 

Cy 

K 

<4 

a§xy 

w  ( kg /nr’ ) 

800 

sin 

15 

l 

3 

0 

sin 

-25 

1.5 

0 

sin 

80 

0.4 

Variable 

hz 

<t>; 

Cz 

faxz 

<4 

axz 

hyz 

<4 

ayz 

P  {kg/m3) 

sin 

-0.12 

1.5 

0 

sin 

0.09 

1.5 

COS 

0.09 

i 

u  (m/s) 

cos 

-18 

0.5 

0 

sin 

60 

0.8 

sin 

60 

0.4 

v  (m/s) 

sin 

-30 

1.25 

0 

cos 

120 

0.7 

sin 

120 

0.8 

w  (m/s) 

cos 

35 

1 

0 

sin 

-80 

0.6 

cos 

-80 

0.8 

p  (. N/m 2) 

cos 

-0.35  x  105 

1 

3 

0 

cos 

0.3  x  105 

1 

cos 

0.3  x  105 

0.5 

Table  6.8:  Additional  Manufactured  solution  coefficients  for  3-D  supersonic  Navier- 
Stokes  cases 
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Figure  6.1:  Grid  nodes  and  control  volumes  on  a  uniform  nodal  mesh 
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Figure  6.2:  Grid  nodes  and  control  volumes  on  a  nodal  mesh  modified  to  achieve  uni¬ 
form  control  volumes 
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Figure  6.3:  Discretization  error  L2  norms  in  the  2-D  supersonic  Euler  case  with  varying 
grid  resolution 


Figure  6.4:  Discretization  error  Z^norm  in  the  2-D  supersonic  Euler  case  with  no  TVD 
limiter  for  varying  grid  resolution 
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Figure  6.5:  Discretization  error  L2  norms  in  the  3-D  supersonic  Euler  case  with  varying 
grid  resolution 


Figure  6.6:  Discretization  error  Lonorm  in  the  2-D  subsonic  Euler  case  with  varying 
grid  resolution 
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Figure  6.7:  Discretization  error  Z^norm  in  the  3-D  subsonic  Euler  case  with  varying 
grid  resolution 


Figure  6.8:  Discretization  error  Li  norm  in  the  2-D  subsonic  Navier-Stokes  case  with 
varying  grid  resolution 
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Figure  6.9:  Discretization  error  Z^norm  in  the  3-D  subsonic  Navier-Stokes  case  with 
varying  grid  resolution 


Figure  6.10:  Discretization  error  Li  norm  in  the  2-D  supersonic  Navier-Stokes  case 
using  original  “first  order”  extrapolation  at  outflow  boundaries  with  varying  grid  reso¬ 
lution 
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Figure  6.11:  Discretization  error  Li  norm  in  the  2-D  supersonic  Navier-Stokes  case 
using  zero^1  order  extrapolation  at  outflow  boundaries  with  varying  grid  resolution 


Figure  6.12:  Discretization  error  Li  norm  in  the  2-D  supersonic  Navier-Stokes  case 
using  corrected  first  order  extrapolation  at  outflow  boundaries  with  varying  grid  reso¬ 
lution 
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Figure  6.13:  Discretization  error  Li  norm  in  the  2-D  supersonic  Navier-Stokes  case 
using  corrected  first  order  extrapolation  at  outflow  boundaries  and  no  TVD  limiter  with 
varying  grid  resolution 


Figure  6.14:  Discretization  error  Li  norm  in  the  2-D  supersonic  Navier-Stokes  case 
using  frozen  outflow  boundaries  with  varying  grid  resolution 
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Figure  6.15:  Discretization  error  Li  norm  in  the  2-D  supersonic  Navier-Stokes  case 
using  frozen  outflow  boundaries  and  the  Koren  limiter  with  varying  grid  resolution 
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Figure  6.16:  Discretization  error  Li  norm  in  the  3-D  supersonic  Navier-Stokes  case 
with  varying  grid  resolution 
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Figure  6.17:  Discretization  error  Linorm  with  varying  grid  resolution  for  the  3-D 
subsonic  Navier-Stokes  case  using  multi-block  grids  with  double-fringe  overlapping 
boundaries 


Figure  6.18:  Discretization  error  Li  norm  with  varying  grid  resolution  for  the  2-D  sub¬ 
sonic  Euler  case  using  a  third  order  scheme 
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Figure  6.19:  Discretization  error  Li  norm  with  varying  grid  resolution  for  the  2-D  sub¬ 
sonic  Euler  case  using  a  third  order  scheme  with  no  limiter 


Figure  6.20:  Discretization  error  Linorm  with  varying  grid  resolution  for  the  2-D  sub¬ 
sonic  Euler  case  using  a  third  order  scheme  with  the  Koren  limiter 
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Figure  6.21:  Discretization  error  Li  norm  with  varying  grid  resolution  for  the  2-D  sub¬ 
sonic  Euler  case  using  a  fourth  order  scheme 


Figure  6.22:  Discretization  error  Linorm  with  varying  grid  resolution  for  the  2-D  sub¬ 
sonic  Euler  case  using  a  fifth  order  scheme 
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Figure  6.23:  Grid  with  density  contours  for  the  subsonic  Euler  test  case  on  a  3-D 
curvilinear  grid 


Figure  6.24:  Discretization  error  Linorm  with  varying  grid  resolution  for  the  3-D  sub¬ 
sonic  Euler  case  using  a  curvilinear  grid 
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Chapter  7 

Significant  Code  Infrastructure 
Improvements 


As  stated  above,  the  ultimate  objective  of  this  project  is  to  obtain  an  improved  capa¬ 
bility  to  perform  time  accurate  simulations  of  large  ground  test  facilities.  Because  of 
the  geometric  complexity  and  large  size  of  these  facilities,  this  requirement  can  only  be 
met  with  a  full-featured  production  CFD  solver  with  “higher”  resolution  features  added 
to  allow  the  transport  of  the  unsteady  features  of  interest  over  the  required  distance.  Be¬ 
cause  of  this,  single-purpose  research  codes,  while  useful  as  a  technology  testbed,  were 
not  sufficient.  Therefore,  transfer  of  any  new  technology  to  a  “production”  CFD  tool 
was  essential  for  the  success  of  the  project.  Much  effort  was  invested  into  building  up 
the  capability  of  a  production  solver  in  order  to  pave  the  way  for  the  desired  large-scale 
unsteady  wind  tunnel  simulations.  As  mentioned  above,  the  Wind-US  code  was  used 
for  this  work. 


7.1  Second  Order  Implicit  Time  Marching 

One  of  the  factors  which  has  limited  Wind-US’s  applicability  in  the  realm  of  unsteady 
flows  has  been  the  implicit  time  marching  procedure.  The  default  algorithm,  which 
came  from  NASTD,  was  a  first  order  backward  Euler  scheme.  This  scheme  is  unsuited 
to  time  accurate  simulations,  and,  indeed,  it  was  never  intended  that  it  be  used  in  such 
a  fashion.  In  NASTD,  unsteady  simulations  always  used  one  of  the  available  Runge- 
Kutta  explicit  time  marching  procedures.  For  the  current  project,  however,  the  ability 
to  run  time-accurate  simulations  of  complex  geometries  with  very  large  time  steps  is 
required,  and  this  necessitates  an  implicit  algorithm. 

To  illustrate  the  problems  with  the  original  implicit  time  marching  algorithm  in 
Wind,  a  simple  inviscid  vortex  convection  case  has  been  run.  This  case  was  obtained 
from  the  suite  of  test  cases  that  are  distributed  with  the  OVERFLOW  solver[40].  A 
vortex  in  a  uniform  mean  flow  is  simulated  on  a  uniform  80x80  Cartesian  grid  (not 
including  boundary  points)  covering  10  ft  on  each  side.  The  boundary  conditions  are 
periodic  so  that  as  the  vortex  is  swept  out  the  back  of  the  domain,  it  reappears  at  the 
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inflow.  The  freestream  Mach  number  was  0.5,  the  total  temperature  was  525°R,  and 
the  total  pressure  was  0.937  psi.  A  non-dimensional  time  step  of  0.1  was  used,  which 
translates  to  a  dimensional  step  of  9.124  x  10  5  seconds  and  200  time  steps  per  flow¬ 
through  time.  The  default  second  order  Roe  spatial  scheme  and  ADI  matrix  inversion 
algorithms  were  employed. 

Fig.  7.1  shows  density  contours  for  both  the  initial  conditions  and  after  five  flow¬ 
through  times.  From  the  figure,  it  is  obvious  that  the  default  scheme  is  inappropriate  for 
unsteady  time  marching.  By  the  end  of  five  flow-through  times,  significant  dissipation 
and  dispersion  errors  are  present. 

As  a  first  step  toward  improving  the  unsteady  implicit  capability,  a  Newton  itera¬ 
tion  procedure  was  added  relatively  early  on  in  the  process  of  developing  Wind.  This 
improved  the  time  marching  ability,  but  persistent  problems  were  present  in  the  im¬ 
plementation  which  have  only  recently  been  resolved  (for  Wind-US).  Fig.  7.2  shows 
the  results  for  the  same  test  case  using  Newton  iterations  in  conjunction  with  the  first 
order  time  marching  algorithm.  For  this  case,  five  sub-iterations  were  used  at  each  time 
level.  The  addition  of  Newton  iterations  reduces  the  dispersion  error  considerably,  but 
the  dissipation  in  the  scheme  is  still  unacceptably  large. 

A  second  order  implicit  time  marching  algorithm  has  now  been  implemented  in  the 
Wind-US  code.  This  technique  is  applicable  either  with  or  without  Newton  iterations. 
At  present,  all  of  the  implicit  schemes  in  Wind,  except  for  the  MacCormack  implicit 
scheme,  can  be  used  with  this  second  order  capability.  Figure  7.3  shows  results  for 
the  same  case  using  the  second  order  scheme  without  any  Newton  iterations.  As  the 
figure  clearly  shows,  the  dissipation  is  greatly  reduced,  but  the  dispersion  error  is  still 
unacceptably  large. 

Finally,  a  run  was  performed  using  the  second  order  algorithm  with  Newton  iter¬ 
ations.  Again,  five  sub-iterations  were  used  at  each  time  level.  Figure  7.4  shows  the 
resulting  density  contours.  In  this  case,  not  only  has  the  dissipation  been  dramatically 
reduced  by  the  second  order  scheme,  compared  to  the  first  order  algorithm,  but  the 
Newton  iteration  procedure  has  dramatically  reduced  the  dispersion  error.  This  result, 
while  still  being  more  dissipative  than  the  high  resolution  schemes  typically  used  for 
the  most  demanding  unsteady  simulations,  is  clearly  a  significant  improvement  over 
the  original  algorithm  in  Wind-US.  Note  that  the  spatial  differencing  used  for  these 
tests  was  the  default  second  order  scheme.  Wind-US  has  the  ability  to  run  schemes  of 
up  to  a  nominal  fifth  order  in  space. 

7.2  Improved  In  o  w  Boundary  Conditions 

In  an  interior  flow  case,  it  is  often  convenient  to  specify  total  conditions  at  an  inflow 
boundary.  In  the  course  of  FY2003,  a  user  of  the  Wind  code  reported  problems  in 
using  the  “hold  totals”  inflow  boundary  condition.  As  a  result  of  this  problem,  the 
inflow  condition  was  rewritten  in  a  greatly  improved  form  that  is  both  more  stable  and 
more  accurate  than  the  previous  version.  While  this  project  did  not  pay  for  that  work, 
it  does  stand  to  benefit  from  it,  and  thus  a  brief  overview  is  presented. 

The  idea  of  both  the  old  and  the  new  boundary  conditions  is  the  same.  The  user 
specifies  a  total  temperature,  total  pressure  and  flow  angle  at  the  boundary  (which  can 
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vary  from  point  to  point).  These  quantities  are  held  fixed  while  the  R2  Riemann  invari¬ 
ant  (u  —  y~j  )  is  extrapolated  from  the  interior.  The  previous  implementation,  how¬ 
ever,  used  a  Newton  iteration  to  simultaneously  solve  for  the  temperature  and  velocity 
at  the  boundary  point.  This  iterative  procedure  would  at  times  not  converge,  resulting 
in  instability  in  the  overall  solution.  The  new  implementation  first  solves  a  quadratic 
equation  for  the  speed  of  sound  based  on  substituting  the  normal  velocity  between  the 
Toand  R2  equations.  This  yields  the  temperature,  and  the  velocity  can  then  be  found 
from  the  definition  of  f?2- 

The  new  approach  is  generally  the  same  as  that  taken  in  the  OVERFLOW  code  [40]. 
In  Wind-US,  however,  the  Riemann  invariant  is  extrapolated  using  a  first  order  algo¬ 
rithm,  rather  than  the  more  usual  zeroth  order.  This  seemingly  simple  change  makes 
a  noticeable  improvement  in  the  total  pressure  distribution  interior  to  the  boundary,  as 
shown  in  Fig.  7.5  for  a  representative  duct  case.  It  is  especially  telling  that  one  cannot 
even  obtain  a  converged  answer  at  all  using  the  old  boundary  conditions  unless  one  sets 
up  the  flow  with  some  other  inflow  boundary  treatment,  whereas  the  new  conditions  are 
stable  enough  to  be  used  from  the  start.  As  is  typical  of  characteristic-based  boundary 
conditions,  there  are  fluctuations  in  the  total  conditions  immediately  downstream  of  the 
boundary.  With  the  new  implementation,  however,  the  values  eventually  return  to  the 
freestream  conditions,  whereas  with  the  previous  boundary  conditions  (and  indeed  the 
new  implementation  as  well,  unless  a  first  order  extrapolation  of  the  invariant  is  used) 
the  flow  never  returns  to  the  freestream  conditions.  This  is  illustrated  in  Fig.  7.6,  which 
shows  the  normalized  total  pressure  along  the  center  line. 

7.3  Compressible  version  of  SST  turbulence  model 

Although  the  majority  of  the  flowfield  in  the  16-S  tunnel  test  case  is  nearly  incompress¬ 
ible,  it  is  anticipated  that,  at  least  in  some  cases  of  facility  unsteadiness,  the  effects  of 
compressibility  could  become  important.  To  address  this  concern,  a  density  weighted 
version  of  the  SST  model  was  written  to  complement  the  original  model,  which  used 
an  incompressible  formulation.  This  version  is  implemented  in  conservative  form  and 
includes  the  compressibility  corrections  of  Suzen  and  Hoffmann[41],  In  addition  to  the 
current  application,  these  terms  are  expected  to  be  especially  helpful  in  the  simulation 
of  free  shear  layers  at  transonic  speeds  (and  higher)  and  other  flows  where  the  effects  of 
compressibility  are  expected  to  be  significant.  The  compressibility  corrections  include 
a  pressure  dilatation  model  and  a  turbulent  Mach  number  based  model  of  compressible 
dissipation. 

To  demonstrate  the  effect  of  the  various  changes  to  the  SST  model,  a  simple  axi- 
symmetric  jet  case  was  performed.  This  case  involves  a  "submerged"  turbulent  su¬ 
personic  jet  emanating  from  an  axi-symmetric  convergent-divergent  Mach  2.2  nozzle. 
This  nozzle  was  first  studied  by  J.  M.  Eggers  in  1962[42],  The  working  fluid  was  air, 
and  the  nozzle  was  operated  at  the  pressure  ratio  corresponding  to  perfect  expansion. 
The  nozzle  plenum  conditions  were  set  to  a  total  temperature  of  525  °R  and  the  total 
pressure  was  162.2psia.  The  ambient  temperature  was  252°R  and  the  ambient  pressure 
was  14.7psia. 

The  simulated  centerline  axial  velocity  is  plotted  in  Fig.  7.7  for  four  different  cases 
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as  well  as  the  experimental  data  of  Eggers.  The  first  two  curves  on  the  plot  (desig¬ 
nated  SST1  and  SST2  on  the  legend)  lie  directly  on  top  of  each  other.  They  are  for 
the  original  version  of  the  SST  model  and  the  new  conservative  version  with  no  com¬ 
pressibility  corrections  added.  As  expected,  they  are  overly  dissipative  of  the  velocity 
downstream  of  the  potential  core.  The  curve  designated  “SST2-CD"  shows  the  effect 
of  the  compressible  dissipation  correction,  while  the  curve  labeled  “SST2-CD-PD” 
shows  the  same  case  with  pressure  dilatation  and  compressible  dissipation  added  in. 
With  the  addition  of  the  compressibility  corrections,  the  velocity  decay  is  much  better 
predicted  aft  of  the  potential  core,  but  the  length  of  the  potential  core  is  increasingly 
over-predicted.  This  is  in  line  with  the  observations  of  other  researchers.  [43]  Thus,  this 
capability,  while  helpful  in  some  cases,  must  be  used  with  caution. 

7.4  Unsteady  Pseudo-Turbulent  In  o  w  BC 

Real  ground  test  facilities  do  not  have  perfectly  uniform  flow.  The  effect  of  unsteady 
inflow  conditions  is  potentially  significant.  For  instance,  Soteriou  and  Ghoniem  [44] 
show  the  development  of  an  incompressible  spatial  mixing  layer  with  and  without  per¬ 
turbation  at  the  inflow  boundary.  Compared  to  the  perturbed  case,  the  development 
of  the  mixing  layer  without  perturbation  is  delayed  significantly  while  the  Kelvin- 
Helmholtz  instability,  activated  by  minute  numerical  errors,  builds  until  the  layer  fi¬ 
nally  rolls  up. 

Researchers  have  taken  several  approaches  to  this  problem.  Sometimes,  as  in  So¬ 
teriou  and  Ghoniem’s  unperturbed  case,  nothing  is  done  to  emulate  turbulence  at  the 
inflow.  For  those  who  wish  to  account  for  inflow  turbulence,  perhaps  the  simplest  ap¬ 
proach  is  to  add  random  “white  noise”  to  the  mean  flow  profiles  at  the  inflow  boundary. 
Comte  et  al.  [45]  did  this  for  a  plane  mixing  layer  and  again  [46]  for  a  plane  wake. 
This  approach  has  significant  drawbacks,  however,  because  this  incoming  “turbulence” 
has  no  correlation  in  space  or  time,  which  is  unphysical. 

Another  approach  is  to  perturb  the  velocity  at  one  or  more  discrete  frequencies. 
These  frequencies  are  usually  based  on  the  most  unstable  frequency  (and  its  subhar¬ 
monics)  as  predicted  by  linear  stability  analysis.  Examples  of  this  approach,  in  addi¬ 
tion  to  the  previously  mentioned  incompressible  work  by  Soteriou  and  Ghoniem[44], 
are  compressible  mixing  layer  computations  by  Sekar  and  Mukunda  [47]  and  Lele  [48]. 
This  approach  is  quite  effective,  especially  if  one  specifically  desires  to  simulate  a  gen¬ 
uinely  forced  flow,  such  as  in  the  experimental  work  by  LeBoeuf  [49].  The  capability 
to  do  this  already  exists  in  the  Wind-US  code.  Real  turbulence,  however,  does  not  gen¬ 
erally  confine  itself  to  a  few  discrete  frequencies;  therefore,  if  a  more  realistic  inflow  is 
desired,  something  more  sophisticated  must  be  done. 

Another  technique  for  introducing  turbulent  fluctuations  at  an  inflow  boundary  is  to 
use  a  separate  temporal  LES  or  DNS,  running  in  parallel  to  the  spatial  case,  to  generate 
the  inflow  turbulence  [50].  Where  it  is  possible  to  use  such  an  approach,  it  would  seem 
to  be  an  ideal  solution.  Unfortunately,  it  is  computationally  very  expensive.  A  variation 
of  this  technique  is  sometimes  used  in  channel  flows  in  which  the  turbulence  from  a 
given  point  downstream  (where  the  flow  is  assumed  to  be  fully  turbulent)  is  rescaled 
and  introduced  at  the  inflow  boundary.  This  method  can  be  quite  successful,  but  is 
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limited  in  its  application  to  cases  where  some  form  of  scale  similarity  law  can  be  used. 
It  would  be  very  difficult,  if  not  impossible  to  apply  such  a  rescaling  method  to  a  fully 
general  case. 

A  final  method  (or  family  of  methods)  for  adding  turbulence  at  the  inflow  boundary 
attempts  to  emulate  the  full  spectrum  of  turbulence  up  to  the  resolution  of  the  grid.  In 
this  method,  measured  or  approximated  turbulent  intensities  at  the  inflow  are  combined 
with  information  on  the  turbulence  energy  spectrum  (usually  an  assumed  functional 
form)  to  derive  time  varying  velocity  perturbations.  Such  an  approach  has  been  used  by 
Oh  and  Loth  [51]  for  supersonic  shear  layers;  Schwer,  Tsuei,  and  Merkle  [52]  for  low 
speed  reacting  mixing  layers;  and  Lee,  Lele,  and  Moin  [53]  for  spatially  developing 
isotropic  turbulence.  More  recently.  Batten,  Goldberg,  and  Chakravarthy  [54]  have 
developed  a  related  technique  and  applied  it  to  a  variety  of  problems. 

The  present  work  will  use  a  variation  of  the  method  of  Lee,  Lele,  and  Moin.  The 
current  implementation  is  designed  to  be  relatively  general,  but  comparatively  inex¬ 
pensive.  The  drawback  is  that  it  lacks  the  sophistication  of  the  more  complex  methods 
(such  as  Batten’s  or  recent  work  by  Keating  and  Piomelli[55]).  The  current  method 
is  best  applied  to  flows  where  there  is  a  strong  initiator  of  turbulent  growth  within  the 
flowfield  itself,  such  that,  while  the  incoming  turbulent  fluctuations  may  be  an  impor¬ 
tant  trigger,  the  fine  details  of  the  incoming  turbulence  are  not  critical.  An  example  of 
such  a  case  is  the  aforementioned  Samimy  and  Elliott  mixing  layer,  where  the  pres¬ 
ence  of  three  dimensionality  in  the  flow  as  it  passes  the  end  of  the  splitter  plate  appears 
sufficient  to  allow  the  development  of  a  turbulent  mixing  layer. 

Conceptually,  the  present  approach  to  inflow  turbulence  may  be  described  as  com¬ 
puting  a  “box”  of  frozen  turbulence  which  is  convected  by  the  mean  flow  into  the  com¬ 
putational  domain  (see  figure  7. 8). The  first  step  in  this  process  is  to  calculate  a  pseudo¬ 
random  divergence-free  isotropic  velocity  field  which  conforms  to  a  prescribed  energy 
spectrum.  This  field  is  computed  initially  in  a  cubic  domain  that  large  enough  to  en¬ 
compass  the  inflow  boundary  and  is  assumed  periodic  in  all  directions.  It  is  computed 
on  a  uniform  Cartesian  grid  (the  default  is  323,  but  can  be  specified  in  the  input  file). 
Given  the  assumed  energy  spectrum,  which  is  discussed  below,  and  a  sufficiently  low 
peak  wavenumber  (based  on  box  size),  this  level  of  resolution  is  sufficient  to  capture 
the  majority  of  the  energy-containing  wavenumbers,  but  it  does  mean  that  in  regions 
where  grid  points  are  highly  concentrated,  the  computational  domain  is  able  to  resolve 
much  higher  wavenumbers  than  are  calculated  in  the  inflow  turbulence. 

The  method  for  computing  this  solenoidal  isotropic  field  is  as  follows.  One  first 
computes  a  stream  function  energy  distribution  in  spectral  space: 


In  the  above  equation,  E(k)  is  the  desired  energy  distribution  of  the  velocity  field.  The 
use  of  the  above  expression  results  in  a  velocity  field  that  approximately  matches  the 
desired  spectrum.  It  is  assumed  that  the  turbulent  energy  is  distributed  evenly  through¬ 
out  each  wavenumber  shell.  Thus  E  (k)  can  be  expressed  as  a  function  of  the  magnitude 
of  the  wavenumber  vector  ( k  =  Ikl).  The  assumed  form  of  the  energy  spectrum  is  the 
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same  as  that  which  has  been  used  previously  for  similar  purposes  by,  among  others, 
Lee,  Lele,  and  Moin  [56]  [53]: 

E  (k)  =  k4exp  ^—2  ^  (7.2) 

For  convenience,  the  spectrum  is  normalized  such  that: 


Enorm  (&)  dk  —  1 

0 

Thus,  for  the  functional  form  given  in  equation  7.2  above, 

‘"“p(-2(i)2) 

Enorm  (k)  —  3  r—  , 

64  v  2  nkp 

The  non-dimensional  wavenumber  at  which  this  function  reaches  its  peak,  kp,  can  be 
specified  by  the  user  either  directly  or  as  a  dimensional  value  (per  foot). 

Next,  the  phase  angles  are  chosen  randomly  (between  -K  and  n).  Using  this  infor¬ 
mation,  a  stream  function, "T,  ,  may  be  computed  in  spectral  space.  Using  a  Fourier 

transform,  this  leads  to  a  physical  stream  function,  'F,  (x),  from  which  velocities  can 
be  computed. 


(7.3) 


(7.4) 


d'f',  3T\, 

dy  dz 

(7.5) 

d'Vx  d'F- 

dz  dx 

(7.6) 

dVy  dx¥x 
dx  dy 

(7.7) 

If  the  field  is  to  be  absolutely  divergence  free,  it  is  necessary  to  perform  this  calcu¬ 
lation  using  the  same  algorithm  for  taking  derivatives  as  the  Navier-Stokes  solver  will 
use.  The  velocity  field  thus  calculated  will  be  truly  divergence-free  only  with  respect 
to  the  derivative  algorithm  that  was  used  to  compute  it.  For  example,  a  field  that  is 
divergence-free  with  respect  to  spectral  derivatives  is  not  solenoidal  with  respect  to  a 
fourth-order  central  difference  derivative  algorithm  or  any  other  non-spectral  deriva¬ 
tive  scheme.  In  the  current  work,  a  fourth  order  central  difference  scheme  is  used  for 
the  derivatives.  Since  the  code  solves  the  compressible  Navier-Stokes  equations,  the 
solenoidal  component  introduced  by  the  turbulent  inflow  due  to  the  use  of  a  different 
numerical  scheme  on  the  interior  should  not  cause  problems. 

Because  the  energy  spectrum  of  the  stream  function  was  only  approximate  (for 
non-spectral  derivatives),  the  velocity  field  must  be  re-scaled  in  order  for  its  energy 
spectrum  to  match  the  target  spectrum  exactly.  This  is  done  by  converting  the  physical 
velocity  field  back  into  spectral  space,  where  its  energy  spectrum  is  calculated.  The 
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spectral  velocity  field  is  then  scaled  so  that  the  energy  in  each  wavenumber  shell  will 
match  the  value  predicted  by  the  desired  energy  spectrum.  The  scaled  field  is  then 
transformed  back  into  physical  space.  As  long  as  the  derivative  algorithm  that  is  used 
is  either  spectral  or  a  uniform  operator  in  space  (in  this  instance,  a  fourth-order  central 
difference  scheme),  the  scaled  field  will  retain  the  solenoidal  properties  of  the  original. 

In  order  to  make  use  of  the  pseudo-turbulent  velocity  field  thus  created,  it  is  as¬ 
sumed  that  the  turbulence  is  convected  into  the  domain  at  a  constant  mean  velocity. 
This  allows  one  to  interpolate  the  velocity  fluctuations  from  the  uniform 
grid  on  which  they  were  generated  onto  the  actual  inflow  plane.  This  field,  however, 
does  not  take  into  account  the  experimentally  determined  distribution  of  turbulent  in¬ 
tensity,  nor  does  it  take  into  account  the  variation  of  resolved  and  unresolved  turbulence 
due  to  differences  in  grid  cell  sizes. 

Since  the  energy  spectrum  is  known  at  the  inflow  plane,  it  is  relatively  simple  to 
compute  the  fraction  of  energy  that  is  resolved  at  a  given  grid  location  (from  the  mesh 
size)  and  modify  the  velocity  perturbations  accordingly.  Also,  if  walls  are  present, 
the  user  can  specify  a  desired  boundary  layer  height,  and  the  velocity  fluctuations  are 
scaled  to  roughly  account  for  their  presence. 

For  supersonic  inflow,  perturbations  can  be  directly  added  to  mean  flow.  For  sub¬ 
sonic  flow,  a  less  direct  method  of  injecting  the  perturbations  must  be  used.  First, 
components  of  velocity  tangential  to  the  boundary  are  directly  specified.  The  normal 
perturbations,  however,  are  introduced  indirectly  by  modifying  the  local  total  temper¬ 
ature  and  total  pressure,  and  using  a  characteristics-based  boundary  condition  to  com¬ 
pute  the  normal  velocity  and  the  thermodynamic  variables. 

In  order  to  prevent  the  introduction  of  spurious  periodicity  to  the  solution,  the  ran¬ 
dom  phase  angles  are  updated  at  irregular  intervals  and  the  velocity  field  is  recom¬ 
puted.  Changing  all  of  the  phase  angles  simultaneously  would  result  in  a  pseudo- 
turbulent  field  which  had  no  correlation  to  the  one  used  just  the  previous  iteration. 
Therefore,  only  the  phase  angles  associated  with  one  wavenumber  shell  are  modified 
at  any  given  time.  Thus,  at  any  given  point  where  the  phase  angles  are  being  modified, 
the  changes  are  relatively  small,  with  a  strong  correlation  to  the  previous  field.  The 
changes  are  spaced  in  such  a  way,  however,  that  in  a  given  flow-through  time  of  the 
pseudo-turbulent  box,  all  of  the  wavenumber  shells  are  thus  “scrambled"  at  least  once, 
resulting  in  a  completely  new  realization. 

Figure  7.9  gives  some  indication  of  the  sort  of  unsteady  field  which  results  from 
the  application  of  this  boundary  condition.  Contours  of  the  tangential  (v)  component  of 
velocity  are  shown  for  supersonic  and  subsonic  2-D  jet  flows.  The  jet  enters  the  domain 
on  the  left  side,  and  covers  the  middle  fifty  percent  of  the  boundary.  The  supersonic 
jet  has  a  mean  core  Mach  number  of  1.209  compared  to  the  surrounding  air,  which  is 
moving  at  Mach  number  1.009.  The  subsonic  jet  case  has  the  freestream  moving  at 
M  =  0.7  and  the  mean  velocity  of  the  jet  at  M  =  0.9.  While  this  case  does  not  properly 
resolve  the  shear  layers  at  the  edge  of  the  jet,  it  clearly  shows  the  effect  of  the  unsteady 
inflow,  as  well  as  the  code’s  ability  to  carry  that  unsteady  flow  downstream. 
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7.5  Other  Fixes  and  Enhancements 

Additional  fixes  and  improvements  to  the  Wind-US  code  have  been  made  during  the 
course  of  this  work.  Although  many,  if  not  most,  of  them  were  not  paid  for  by  AFOSR, 
the  Facility  Unsteadiness  Research  project  benefits  from  the  improved  computational 
capability.  Some  of  the  more  relevant  fixes  are  now  noted  in  brief. 

Wind-US  (and  Wind  before  it)  had,  at  least  nominally,  simple  steady-state  screen 
models.  Unfortunately,  the  implementation  of  these  models  had  serious  bugs  that  pre¬ 
vented  them  from  functioning  properly.  These  problems  have  now  been  fixed  and  the 
model  appears  to  work  as  advertised.  This  provides  an  important  starting  point  for  any, 
more  advanced,  models  which  may  be  developed  in  future. 

Persistent,  but  intermittent,  problems  have  been  encountered  in  the  Wind  code  for 
quite  some  time  in  certain  areas  that  were  critical  to  the  successful  execution  of  this 
project.  Specifically,  the  Newton  iteration  and  block-to-block  coupling  algorithms  had 
numerous  problems  related  to  memory  addressing,  parallel  processing,  and  unantici¬ 
pated  special  case  situations.  At  irregular  intervals,  these  would  cause  the  code  to  fail 
without  warning  (and  often  without  any  error  message  at  all).  Thanks  to  a  significant 
rewrite  which  allows  Wind,  and  now  Wind-US,  to  be  run  with  very  restrictive  error 
checking  under  a  debugger,  these  algorithms  have  been  stabilized  (and  many  other 
problems  fixed  as  well). 

For  some  time,  there  has  been  a  desire  for  the  capability  to  run  Wind-US  with  dou¬ 
ble  precision  variables.  A  variety  of  applications,  including  the  Facility  Unsteadiness 
Research  project,  were  seen  as  possible  beneficiaries  of  such  a  capability,  but  the  work 
involved  in  rewriting  the  code  was  seen  as  too  great,  and  the  task  was  never  undertaken 
until  recently,  when  an  overall  rewrite  of  the  memory  allocation  to  use  FORTRAN 
90/95  constructs  was  performed  (mostly  by  Boeing).  At  that  time,  a  capability  was 
introduced  to  specify  the  precision  of  the  variables  at  compile  time.  Although  this 
capability  is  fully  functional,  work  is  still  ongoing  to  ensure  that  the  precision  of  the 
solution  is  not  corrupted  by  lower  precision  constants  or  variables  being  used  inadver¬ 
tently  during  the  solution  process. 

One  area  which  has  seen  an  immediate  impact  of  the  new  capability  is  the  Method 
of  Manufactured  Solutions  verification  process;  it  is  now  possible  to  converge  resid¬ 
uals  at  least  a  further  six  orders  of  magnitude  beyond  what  the  single  precision  code 
can  accomplish.  This  removes  a  significant  unknown  from  the  MMS  process,  because, 
especially  for  the  higher  order  schemes,  the  error  was  often  small  enough  that  the  pre¬ 
cision  of  the  code  prevented  a  fully  converged  realization  of  the  problem.  An  example 
of  this  is  shown  in  Figure  7.10.  The  first  plot  shows  the  error  in  the  total  energy  vari¬ 
able  after  the  solution  has  been  converged  as  much  as  possible  using  a  single  precision 
code.  The  second  plot  is  the  same  case  run  with  the  double  precision  version.  With 
the  improved  precision,  although  the  overall  error  levels  remain  similar,  the  solution 
is  much  smoother,  indicating  that  what  is  being  observed  is  truly  the  error  in  the  nu¬ 
merical  scheme,  rather  than  the  results  of  precision  limitations.  Similar  benefits  are 
expected  to  accrue  for  long  unsteady  runs,  where  tens  and  hundreds  of  thousands  of 
realizations  are  to  be  averaged  in  order  to  obtain  the  desired  statistics.  Without  double 
precision  variables,  such  calculations  are  subject  to  significant  error. 

Another  area  of  work  has  been  the  verification  and  validation  of  the  various  ca- 
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pabilities  implemented  in  Wind-US.  Some  of  this  work  is  documented  elsewhere  (in 
Chapter  6).  A  related  issue  that  has  been  a  persistent  problem  for  NPARC  Alliance 
developers  has  been  the  frequency  with  which  a  fix  for  one  bug  introduces  a  problem 
in  another  part  of  the  code  which  heretofore  worked.  To  address  this  problem,  an  effort 
was  initiated  to  begin  gathering,  organizing,  and  computing  baseline  solutions  for  a  set 
of  test  cases  that  would  exercise  as  many  features  of  the  code  as  could  be  managed. 

The  ultimate  goal  is  to  have  suites  of  test  cases  which  would  cover  the  full  range  of 
physics,  models,  and  code  options  available  in  the  code.  Along  with  the  cases,  scripts 
were  written  to  automatically  traverse  a  directory  tree,  run  each  case  that  is  encoun¬ 
tered,  and  compare  with  the  baseline  solution.  Currently,  cases  have  been  identified 
which  span  a  significant  portion  of  the  structured  solver’s  capabilities.  In  addition, 
some  unstructured  cases  have  been  identified,  but  the  coverage  is  less  complete  in  this 
area. 

The  test  cases  are  divided  into  several  “suites”  of  problems.  The  most  commonly 
used  one  is  a  set  of  problems  that  run  relatively  quickly  (in  serial  mode).  The  purpose 
of  this  suite  is  to  test  the  basic  mechanics  of  the  code.  A  second  test  suite  consists 
of  cases  which  run  in  parallel  in  order  to  ensure  that  the  algorithms  still  function  in 
that  mode.  A  third  suite  consists  of  cases  which  are  much  larger,  and  they  test  the 
capability  of  the  code  to  handle  complex  file  and  memory  requirements,  as  well  as 
cases  with  many  options  interacting  with  each  other.  The  above  tests  generally  run 
only  for  a  few  iterations;  another  suite  includes  cases  which  are  run  to  convergence. 

This  ensures  that  the  accuracy  of  the  code  results  have  not  been  compromised. 

Note  that,  although  there  is  some  overlap,  these  suites  are  distinct  from  the  “NPARC 
Alliance  CFD  Verification  and  Validation  Website”  (http://www.grc.nasa.gov/WWW/wind/valid/validation.html). 
The  web  site  contains  many  good  test  cases,  along  with  results.  The  Validation  Web¬ 
site,  however,  treats  each  case  as  an  individual  exercise;  there  is  no  provision  made 
for  running  groups  of  cases  en  masse.  The  “test  suites”  are  an  attempt  to  remedy  that 
situation. 

Other  minor  fixes  and  enhancements  have  focused  on  improving  the  ability  to  col¬ 
lect  and  process  the  massive  quantities  of  data  that  a  large  scale  unsteady  simulation 
will  generate.  The  end  result  of  these  enhancements  is  a  code  which  can  now  (at  last) 
begin  to  investigate  unsteady  flow  in  large,  complex  facilities. 

7.6  WICS  Bay  Validation  Case 

The  WICS  (Weapons  Internal  Carriage  and  Separation)  L/D  =  4.5  weapons  bay  [31] 
for  M  =  0.95  and  Re  =  2.5  x  10 6 /ft  has  been  used  as  a  check  case  to  assist  in  the 
validation  of  many  of  the  aforementioned  improvements  to  Wind-US.  This  is  the  same 
case  discussed  earlier  in  Chapter  5,  where  it  was  run  using  the  NXAIR  code.  The  wind 
tunnel  data  were  obtained  in  the  AEDC  four-foot  transonic  wind  tunnel.  The  weapons 
bay  was  18  in.  long,  4  in.  wide,  and  4  in.  deep.  The  computational  geometry  was  a 
flat  plate  that  extended  15  in.  upstream  of  the  bay  to  match  the  experiment  and  25  in. 
downstream  of  the  bay.  The  sides  of  the  computational  grid  extended  50  in.  on  either 
side  of  the  bay  centerline.  Wall  functions  were  employed  for  all  viscous  walls.  The  wall 
spacing  was  chosen  as  0.0075  in.,  which  corresponds  to  a  >’+  of  50  on  the  upstream 
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plate.  The  grid  system  had  1,132,860  points  broken  into  eight  grids.  A  constant  time 
step  of  7.88  x  1 0  6  seconds  was  used  in  the  calculation.  These  cases  were  run  using 
the  second  order  implicit  scheme  with  Newton  iterations  discussed  above.  Results  have 
been  obtained  for  both  the  DES  model  and  the  hybrid  turbulence  model  of  Nichols  and 
Nelson.  The  latter  model  used  the  compressible  form  of  the  Menter  SST  turbulence 
model  discussed  above,  including  the  compressibility  corrections. 

The  calculation  was  run  for  3,000  steps  to  remove  the  initial  transients  and  ensure 
that  the  flow  was  fully  developed  (displaying  significant  three-dimensionality).  Data 
was  then  gathered  for  at  least  8192  time  steps  (the  final  8192  time  steps  being  used  for 
spectral  analysis)  similar  to  the  practice  used  when  running  this  case  with  NXAIR. 

Figure  7.11  shows  Mach  number  contours  from  simulations  using  the  DES  and 
Nichols  and  Nelson  hybrid  models  (at  different  times).  These  results  show  similar  be¬ 
havior  as  those  reported  in  Chapter  5  (Figure  5.16).  Although  the  flowfield  snapshots 
are  taken  from  slightly  different  times  in  the  simulation,  they  are  qualitatively  very 
similar.  Similar  levels  of  vorticity  are  observed  in  both  simulations  as  well,  as  shown 
in  Figure  7.12  (for  the  Z-component).  As  with  the  earlier  NXAIR  runs,  the  turbulent 
viscosity  of  the  DES  model  is  noticeably  lower  than  that  of  the  hybrid  model  (as  shown 
in  Figure  7.13).  This  reflects  the  fact  that,  away  from  walls,  the  DES  model  switches  to 
a  simple  LES  mode  regardless  of  flow  conditions,  whereas  the  hybrid  model  attempts 
to  estimate  the  ratio  of  resolved  to  unresolved  turbulent  kinetic  energy  at  each  point  and 
bases  the  switch  to  LES  mode  on  that.  As  a  result,  the  hybrid  model  is  detecting  re¬ 
gions  away  from  walls  where  the  turbulent  scales  are  small  enough  that  the  grid  cannot 
resolve  a  substantial  fraction  of  the  turbulent  kinetic  energy.  To  compensate,  the  tur¬ 
bulent  viscosity  is  increased.  Despite  this  difference,  both  models  predict  remarkably 
similar  flow  structures.  These  similarities  are  not  confined  to  spatial  quantities.  The 
sound  pressure  level  spectra  at  two  different  points  along  the  centerline  of  the  cavity 
wall  are  plotted  in  Figure  7.14.  Although  some  of  the  details  differ,  the  overall  lev¬ 
els  and  trends  are  virtually  identical.  Visually  comparing  these  plots  with  the  results 
of  the  NXAIR  runs  presented  earlier  indicates  that  both  codes  are  predicting  similar 
structures.  The  sound  pressure  level  spectrum  predicted  by  Wind-US  would  appear  to 
have  a  somewhat  larger  magnitude  than  that  predicted  by  NXAIR  or  observed  in  the 
experiment.  This  may  reflect  a  genuine  difference  between  the  two  codes,  but  it  could 
equally  well  be  due  to  differences  in  post-processing.  Unfortunately,  a  direct  compari¬ 
son  with  the  results  from  NXAIR  and  the  experiment  was  not  possible,  as  the  original 
data  was  not  available  to  the  author. 

Thus,  once  again,  all  else  being  equal,  both  the  DES  and  Nichols  and  Nelson  ap¬ 
proaches  to  hybrid  RANS/LES  turbulence  modeling  yield  functionally  equivalent  re¬ 
sults.  This  is  despite  fairly  significant  differences  in  the  predicted  levels  of  turbulent 
viscosity.  As  before,  it  appears  that  for  cases  such  as  this,  other  factors,  such  as  the 
time  marching  scheme,  may  be  more  important  than  the  particular  hybrid  model  used. 

7.12 
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Figure  7.1:  Initial  and  final  density  contours  from  an  inviscid  vortex  convection  simu¬ 
lation  using  first  order  time  marching 


Figure  7.2:  Initial  and  final  density  contours  from  an  inviscid  vortex  convection  simu¬ 
lation  using  first  order  time  marching  with  Newton  iteration 
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Figure  7.3:  Initial  and  final  density  contours  from  an  inviscid  vortex  convection  simu¬ 
lation  using  second  order  time  marching 


Figure  7.4:  Initial  and  final  density  contours  from  an  inviscid  vortex  convection  simu¬ 
lation  using  second  order  time  marching  with  Newton  iterations 
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(a)  Original  Wind  Infbw  Boundary  Condition 


(b)  Improved  Wind-US  Inlbw  Boundary  Condition 

100 

Figure  7.5:  Normalized  total  pressure  field  variations  with  inflow  boundary  condition 
treatment 
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Figure  7.6:  Effect  of  inflow  boundary  condition  treatment  on  normalized  total  pressure 


Figure  7.7:  Effect  of  compressibility  corrections  on  predictions  of  centerline  velocity 
in  an  axi-symmetric  nozzle 
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Figure  7.8:  Schematic  of  the  pseudo-turbulent  unsteady  inflow  boundary  condition 


(a)  Supersonic  Jet 


(b)  Subsonic  jet 


Figure  7.9:  V  component  of  velocity  in  a  2-D  turbulent  jet  case  with  pseudo-turbulent 
inflow  (flow  from  left  to  right) 
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(a)  Single  Precision  Solver 


(b)  Double  precision  solver 


Figure  7.10:  Effect  of  precision  on  the  error  in  total  energy  in  a  Method  of  Manufac¬ 
tured  Solutions  case  on  a  643  grid 
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(a)  Nichols  and  Nelson  hybrid  model 


Figure  7.11:  Contours  of  Mach  number  in  unsteady  WICS  bay  simulations 
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(a)  Nichols  and  Nelson  hybrid  model 
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(b)  DES  model 


Figure  7.12:  Contours  of  the  Z-component  of  vorticity  in  unsteady  WICS  bay  simula¬ 
tions 


105 


AEDC-TR-05-1 


(a)  Nichols  and  Nelson  hybrid  model 


(b)  DES  model 


Figure  7.13:  Contours  of  turbulent  viscosity  in  unsteady  WICS  bay  simulations 
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(a)  Bottom  wall  X/L  =  0.985 


(b)  Back  Wall  Z/L  =  0.85 


Figure  7.14:  Sound  pressure  level  along  bay  centerline  in  unsteady  WICS  bay  simula¬ 
tions 
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Chapter  8 

Prediction  of  Diffuser  Instability 


It  is  believed  that  the  major  contributing  factor  to  the  unsteadiness  experienced  in  the 
16-S  wind  tunnel  at  AEDC  is  intermittent  instability  in  the  diffuser  section  between 
corner  three  and  corner  four  (see  Fig.  8.1).  This  sort  of  instability  has  been  seen  in  the 
past  in  other  facilities.  It  thus  represents  a  class  of  facility  problems  that  would  have 
to  be  numerically  resolvable  in  order  to  have  a  reasonably  general  facility  unsteadiness 
simulation  capability.  The  following  cases  have  been  undertaken  in  order  to  pave  the 
way  toward  (the  capability  for)  simulations  of  the  full  16-S  tunnel  benchmark  case. 

8.1  Steady  State  Runs 

Experimentally,  flow  in  diffusers  has  been  observed  to  be  completely  attached,  inter¬ 
mittently  separating,  or  continually  separated,  depending  on  the  length  of  the  diffuser 
and  the  spreading  angle.  The  intention  in  this  project  is  to  simulate  various  data  points 
(shown  in  Fig.  8.2)  to  determine  if  CFD  can  match  the  experimentally  observed  behav¬ 
ior.  As  a  first  step  in  this  direction,  several  two-dimensional  axi-symmetric  steady-state 
cases  have  been  run.  The  intention  is  to  determine  whether  or  not  the  Wind-US  code 
has  the  capability  to  capture  the  most  basic  aspects  of  the  physics.  If  the  code  cannot 
capture  at  least  the  steady-state  aspects  of  these  diffusers,  then  it  is  doubtful  that  an 
unsteady  simulation  could  yield  much  useful  information. 

The  results  from  steady  state  runs  were  extremely  encouraging.  Cases  four  and 
five,  which  should  be  in  the  fully  attached  region,  do  in  fact  show  fully  attached  flow 
as  shown  in  Figs.  8.3  and  8.4.  The  cases  which  should  have  large  separation  regions, 
are  also  correctly  predicted  (see  Figs.  8.5  and  8.6).  Most  interesting  of  all  is  the  case 
that  straddles  the  stall  line;  this  case,  when  examined  in  detail  (see  Fig.  8.7),  shows 
incipient  separation,  just  as  it  should.  Thus,  Wind-US  has  been  validated  to  at  least  be 
able  to  capture  the  basic  steady-state  flow  features  in  these  diffusers. 
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8.2  2-D  Simpli  ed  16-S  Tunnel  Model 

In  another  effort  to  assess  the  ability  of  Wind-US  to  meet  the  requirements  for  this 
project,  a  simplified,  2-D  grid  of  the  16-S  tunnel  was  created.  This  grid  runs  from 
“valve  80”  through  the  test  section  (see  Fig.  8.1),  but  is  purely  two  dimensional  and 
does  not  include  any  turning  vanes  or  screens.  Obviously,  in  such  a  configurations, 
one  expects  separation  at  the  corners,  and  the  results  will  bear  little  relation  to  what 
occurs  in  the  actual  tunnel.  The  objective  of  this  case,  however,  is  not  to  model  16-S, 
per  se,  but  to  determine  if  the  code  can  support  unsteadiness  sufficiently  that  the  fluc¬ 
tuations  due  to  the  separation  at  corner  four  is  detectable  in  the  test  section.  The  results 
using  the  default  (second  order  in  space  )  scheme  indicate  that  some  unsteadiness  is 
being  predicted,  and  that  some  of  this  unsteadiness  appears  to  be  surviving  into  the 
test  section.  The  magnitude  of  pressure  variation  with  time  at  points  along  the  wall  of 
the  test  section  is  on  the  order  of  0.03  psi  (out  of  approximately  1.5  psi).  Figure  8.8 
shows  contours  of  the  Z-component  of  vorticity  for  a  case  run  with  the  DES  model 
compared  with  a  steady-state  case  run  with  Menter’s  SST  model.  The  unsteady  case 
clearly  shows  some  vortices  being  shed  at  Corner  3,  but  these  seem  to  be  damped  out 
rather  quickly. 

8.3  Axi-Symmetric  Simpli  ed  16-S  Tunnel  Model 

A  second  simplified  version  of  the  16-S  tunnel  was  created  in  an  effort  to  address 
the  shortcomings  of  the  above  test  case.  The  new  version  uses  a  2-D  axi-symmetric 
grid  which  consists  of  the  diffuser  leg,  stilling  chamber,  and  test  section  attached  in 
a  straight  line.  This  is  roughly  equivalent  to  assuming  that  the  turning  vanes  perform 
perfectly  and  turn  the  flow  with  no  losses  whatsoever.  It  also  changes  the  nozzle  and 
test  section  cross-sections  from  rectangular  to  circular.  One  immediate  result  of  this 
is  that  the  area  ratios  of  the  nozzle  are  no  longer  appropriate  for  Mach  1.6  flow  and 
instead  yield  Mach  numbers  of  roughly  1 .9  in  the  test  section. 

Despite  these  drawbacks,  if,  as  is  believed,  the  diffuser  is  unstable,  this  case  should 
be  sufficient  to  capture  some  degree  of  unsteadiness  and  propagate  it  downstream  to  the 
test  section  (subject  to  the  limitations  of  the  numerics).  Using  this  case,  therefore,  it 
should  be  possible  to  obtain  an  estimate  of  what  a  comparable  3-D  case  would  require 
in  terms  of  computational  resources  (due  to  grid  resolution  requirements  and  time  step 
limitations). 

Mach  number  contours  for  this  case  are  shown  in  Figure  8.9.  From  this  plot,  it 
would  seem  that  no  noticeable  unsteadiness  is  being  predicted  in  this  case.  When 
the  details  of  the  flow  in  the  diffuser  section  are  examined,  however,  one  can  see, 
intermittently,  separation  bubbles  forming,  bursting,  and  then  reforming.  The  intensity 
of  the  separation  is  not  as  extreme  as  that  observed  in  the  subscale  tunnel  experiment, 
but  it  is  very  encouraging  nonetheless.  The  region  of  reverse  flow  during  one  of  these 
episodes  of  separation  is  shown  in  Figure  8.10. 

The  contours  of  the  Z-component  of  vorticity  in  Figure  8.11  show  how  the  pseudo- 
turbulent  fluctuations  being  injected  at  the  inflow  are  propagating  through  the  domain. 
Also  visible  are  the  larger  vortical  structures  being  shed  from  the  diffuser  wall  and 
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propagated  downstream.  The  scale  of  plot  has  been  greatly  truncated  so  as  to  high¬ 
light  the  freestream  turbulence,  which  is  relatively  weak:  inflow  streamwise  turbulent 
intensity  is  two  percent  of  the  freestream  velocity.  Also  note  that  the  behavior  is  dif¬ 
ferent  than  observed  in  Figure  7.9  in  that  long  streaks  of  vorticity  can  be  seen  in  the 
streamwise  direction.  At  least  some  of  these  differences  are  due  to  the  pseudo-turbulent 
inflow  being  designed  for  either  true  2-D  or  full  3-D  simulations,  not  axi-symmetric 
flows  such  as  this  one.  The  varying  circumference  as  one  changes  radial  location  will 
distort  the  development  of  the  unsteady  flow  relative  to  what  one  would  expect  in  a 
pure  two-dimensional  case. 

The  effect  of  these  fluctuations  has  been  measured  at  the  end  of  the  nozzle  (start  of 
the  test  section).  Traces  of  streamwise  velocity  and  static  pressure  at  the  centerline  are 
shown  in  Figures  8.12  and  8.13,  respectively.  Since  the  overall  fluctuations  observed 
in  the  diffuser  were  smaller  and  less  intense  than  that  observed  in  the  tunnel,  it  is  no 
surprise  that  the  magnitude  of  the  fluctuations  observed  in  the  test  section  should  be 
small  as  well.  Nevertheless,  it  is  highly  encouraging  to  see  that  these  relatively  weak 
structures  are  in  fact  propagated  through  the  stilling  chamber  and  can  be  detected  in 
the  test  section  some  two  hundred  feet  downstream. 


Figure  8.1:  Schematic  of  the  AEDC  16-S  Wind  Tunnel 
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Figure  8.2:  Observed  Stability  Characteristics  of  Diffusers 
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Figure  8.3:  Contours  of  velocity  magnitude  predicted  for  the  Case  4  diffuser 
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Figure  8.4:  Contours  of  velocity  magnitude  predicted  for  the  Case  5  diffuser 
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Figure  8.5:  Contours  of  velocity  magnitude  predicted  for  the  Case  7  diffuser 
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Figure  8.6:  Contours  of  velocity  magnitude  predicted  for  the  Case  9  diffuser 
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(a)  Full  domain 


Figure  8.7:  Contours  of  velocity  magnitude  predicted  for  the  Case  6  diffuser 
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(a)  Steady  state  run  (SST  model) 


(b)  Unsteady  run  (DES  model) 


Figure  8.8:  Contours  of  Z  vorticity  in  a  simplified  2-D  model  of  16-S 
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Figure  8.9:  Contours  of  Mach  number  in  a  simplified  axi-symmetric  model  of  16-S 


Figure  8.10:  Regions  of  reverse  flow  (blue)  at  a  given  instant  in  a  simplified  axi- 
symmetric  model  of  16-S 


Figure  8.11:  Contours  of  Z-vorticity  in  a  simplified  axi-symmetric  model  of  16-S 
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Figure  8.12:  Trace  of  streamwise  velocity  in  the  test  section  of  a  simplified  axi- 
symmetric  model  of  16-S 
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Figure  8.13:  Trace  of  static  pressure  in  the  test  section  of  a  simplified  axi-symmetric 
model  of  16-S 
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Chapter  9 

Conclusions  and 
Recommendations 


The  Facility  Unsteadiness  Research  project  has  resulted  in  a  much-improved  ability  to 
use  CFD  to  predict  unsteadiness  in  large  scale  ground  test  facilities.  The  capabilities 
thus  developed  have  been  transitioned  to  the  Wind-US  code  for  wider  use  by  the  Air 
Force  and  other  DoD-related  organizations.  This  government  code  continues  to  be  de¬ 
veloped  and  maintained,  providing  some  assurance  that  these  capabilities  will  continue 
to  be  available  to  the  DoD  community. 

This  work  has  shown,  among  other  things,  that  the  basic  ideas  of  hybrid  RANS/LES 
turbulence  models  are  at  least  in  some  sense  correct.  Specifically,  by  transitioning  from 
a  RANS  modeling  approach  in  regions  that  are  under-resolved  by  the  computational 
grid  to  an  LES-type  model  where  the  grid  is  suitably  dense,  dramatic  improvements 
are  possible  in  the  accuracy  of  unsteady  simulations.  This  is  true  for  the  two  models 
discussed  in  this  work  (DES  and  the  hybrid  model),  but  there  is  not  reason  to  believe  it 
would  not  also  be  true  for  other  approaches.  Furthermore,  it  appears  that,  at  least  for  the 
mixing  layer  simulations  discussed  above,  the  differences  between  the  two  approaches 
were  minimal,  with  no  clear  “winner”  yet  in  sight.  Combined  with  the  analysis  which 
showed  that  boundary  location  and  statistical  gathering  methods  have  as  much  or  more 
effect,  it  becomes  obvious  that  expending  more  time  and  effort  on  tweaking  hybrid 
models  would  be  counterproductive,  at  least  for  the  purposes  of  pursuing  large  scale 
unsteady  simulations.  The  one  exception  to  this  would  be  in  the  area  of  the  transition 
between  RANS  and  LES  modes.  The  difficulties  encountered  with  the  DES  model 
when  a  dense  computational  grid  is  present,  but  little  or  no  unsteadiness  is  encountered 
indicate  that  perhaps  some  modification  to  the  switch  between  RANS  and  LES  modes 
would  be  in  order.  The  hybrid  model  of  Nichols  and  Nelson  has  shown  a  somewhat 
more  robust  response  to  the  same  situation,  but  more  work  is  required  in  order  to  de¬ 
termine  if  that  was  merely  fortuitous  or  if  it  represents  an  advantage  of  this  approach 
over  that  taken  in  the  DES  model. 

The  infrastructure  improvements  documented  above  that  have  taken  place  within 
the  Wind-US  code  have  resulted  in  a  much  more  robust  code  and  a  better  platform  for 
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further  development.  The  improved  time  marching,  in  particular,  has  resulted  in  at  least 
an  order  of  magnitude  improvement  in  turn-around  time  for  unsteady  cases.  The  use 
of  the  Method  of  Manufactured  Solutions  and  the  various  test  suites  has  proved  to  be 
an  extremely  useful  technique  for  maintaining  confidence  in  the  code  as  changes  are 
made. 

The  diffuser  test  cases  that  were  run  have  shown  that  the  Wind-US  code  is  capable 
of  matching  the  observed  separation  in  diffusers  in  a  steady-state  sense.  The  improved 
unsteady  capability  has  also  been  put  to  use  in  simulations  of  the  WICS  bay  and  sim¬ 
plified  versions  of  the  16-S  wind  tunnel  geometry.  The  results  from  the  WICS  bay 
simulations  are  in  reasonable  agreement  with  the  findings  of  other  researchers. 

The  simplified  16-S  simulations  indicate  that  the  instability  can  be  captured  in  such 
a  large  facility  and  propagated  for  significant  distances  while  remaining  detectable. 
This  suggests  that  with  a  full  three-dimensional  grid  with  somewhat  more  resolution  in 
the  center  of  the  duct,  the  associated  3-D  instability  modes  and  3-D  pseudo-turbulent 
inflow  fluctuations  could  resolve  the  large  scale  unsteadiness  in  facilities  such  as  the 
16-S  tunnel.  Using  the  technology  developed  and/or  demonstrated  in  this  work,  a  sim¬ 
ulation  of  the  diffuser,  stilling  chamber,  nozzle,  and  test  section  would  require  a  grid 
with  an  estimated  20-30  million  grid  points.  Running  such  a  large  grid  for  the  required 
105  timesteps  (in  order  to  wash  out  initial  conditions  and  gather  sufficient  data  to  allow 
statistically  meaningful  time  averaging)  would  obviously  require  significant  computa¬ 
tional  resources.  Such  resources  exist,  however,  at  the  DoD’s  Major  Shared  Resource 
Centers,  and  similar  sized  cases  are  not  unheard  of.  Thus,  while  time-accurate  simula¬ 
tions  of  large  facility  unsteadiness  are  not  yet  within  the  reach  of  everyday  production 
CFD,  there  is  now  evidence  that  suggests  these  simulations  are  within  the  reach  of  large 
computer  facilities. 

With  this  demonstrated  basic  unsteady  capability  in  place,  future  work  can  now 
proceed  toward  refining  it  and  implementing  more  specialized  technologies  needed 
for  particular  projects.  For  large  scale  test  facility  simulations  in  particular,  the  addi¬ 
tion  of  verified  high-resolution  numerical  schemes  (better  than  the  current  fifth  order 
upwind-biased  schemes  already  present)  would  be  greatly  beneficial,  as  it  would  result 
in  anywhere  from  two  to  eightfold  improvement  in  grid  resolution  requirements  for 
3-D  cases  and  thus  substantially  reduce  the  cost  of  such  simulations. 
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